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I.  INRTODCUCTION 

The  postdoctoral  fellowship  grant  was  awarded  to  the  principal  investigator  (PI)  for  the  period  of  April 
1,  2003 — ^March  31,  2005.  The  purpose  of  this  investigation  is  to  introduce  a  framework  for  including 
model  parameter  uncertainties  into  prostate  Intensity  Modulation  Radiation  Therapy  (IMRT)  dose 
optimization  so  that  biological  model-based  objective  function  can  be  used  with  improved  confidence  level. 
The  specific  aims  of  the  proposal  are:  (1)  to  establish  a  mathematical  formalism  to  incorporate  model 
parameter  uncertainty  into  IMRT  optimization;  (2)  to  identify  the  clinically  relevant  biological  model 
parameter  variance  range;  and  (3)  to  study  the  prostate  cancer  treatment  planning  including  the  model 
uncertainty  information.  Under  the  generous  support  from  the  U.S.  Army  Medical  Research  and  Materiel 
Command  (AMRMC),  the  PI  has  contributed  significantly  to  the  radiation  treatment  of  prostate  cancer. 
Several  conference  abstracts  and  refereed  papers  have  been  resulted  from  the  support.  The  project  is 
integrated  with  ongoing  radiation  treatment  so  that  the  PI  has  obtained  clinical  experience  while 
accomplishing  the  proposed  projects.  The  preliminary  data  and  clinical  training  obtained  under  the  support 
of  this  grant  has  enabled  the  PI  to  start  as  an  assistant  professor  in  the  Department  of  Radiation  Oncology  at 
a  prestigious  university. 


II.RESEARCH  AND  ACCOMPLISHMENTS 


Adenocarcinoma  of  the  prostate  has  become  the  most  common  malignancy  in  men  in  the  western 
countries.  Options  for  active  management  of  organ-confined  prostate  cancer  include  radical  prostatectomy 
and  definitive  radiotherapy  with  either  external  beams  or  interstitial  brachytherapy.  Intensity  Modulated 
Radiation  Therapy  (IMRT)  is  quickly  replacing  conventional  techniques  for  the  treatment  of  prostate 
cancer.  Most  IMRT  optimization  systems  at  present  use  dose  an^or  dose  volume-based  objective 
functions  S  which  guide  the  IMRT  planning  by  imposing  a  penalty  according  to  the  difference  between  the 
computed  and  prescribed  doses.  A  well-known  drawback  of  the  dose-based  inverse  planning  is  that  the 
nonlinear  dose  response  of  tumor  or  normal  structures  is  not  fully  considered.  A  number  of  mathematical 
models  have  been  developed  over  the  years  to  better  describe  the  biological  effect  of  radiation  and 
considerable  works  have  also  been  done  to  use  these  biological  models  to  construct  more  meaningful 
objective  functions  for  therapeutic  dose  optimization  Generally  speaking,  radiobiological  formalism 
involves  the  use  of  model  parameters  that  are  of  considerable  uncertainty  For  instance,  the 
radiosensitivity  a  of  Webb’s  TCP  model  varies  from  0.157  Gy'^  to  0.090  Gy'Vwhen  model  parameters  were 
fit  to  103  patients’  data  Biological  ‘margins’  have  been  used  to  account  for  the  variability  in  radiation 
sensitivity.  Similar  to  the  use  of  a  safety  margin  to  account  for  the  potential  uncertainties  in  targeting  a 
tumor,  this  method  assigns  more  conservative  radiosensitivity  values  to  the  tumor  or  sensitive  structures  to 
deal  with  the  potential  uncertainty  of  the  parameter^. 

The  purpose  of  this  project  is  to  develop  a  framework  to  include  any  types  of  model  parameter 
uncertainties  to  the  dose  optimization.  In  our  approach,  the  uncertainty  of  a  model  parameter  is  quantified 
by  a  probability  density  function  and  its  influence  is  then  incorporated  into  inverse  planning  through  the  use 
of  a  statistical  inference  theorem Our  model  is  formulated  on  the  equivalent  uniform  dose  (EUD)  ^ : 


EUD  = 


jr'Eo; 


(1) 


for  both  tumor  and  normal  tissues,  where  N  is  the  number  of  voxels  in  the  structure,  Di  is  the  dose 
delivered  to  the  ith  voxel,  a  is  the  tumor  or  normal  tissue-specific  parameter  that  describes  the  dose-volume 
effect.  The  corresponding  objective  function  to  measure  the  goodness  of  a  dose  distribution  is  given  by  ® 


'  .  (2) 

where  the  component  subcore^  may  be  either 
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for  tumors,  or 


for  normal  tissues  and  organs  at  risk  (OARs).  EUDo  is  the  desired  dose  parameter  for  the  target  volume  and 
the  maximum  tolerable  xmiform  dose  for  normal  structures.  Parameter  n  is  akin  to  the  structure  specific 
importance  factor  in  the  conventional  inverse  plaiming  formalism  that  parameterizes  our  tradeoff 
strategy  of  different  structure. 

We  assume  that  Ok  in  the  EUD  model  varies  according  to  a  simple  Gaussian  distribution 
Pni.ak)  =  P'^  exp{-r„[at  -aj'} ,  (5) 

where  ao  is  the  mean  value,  P\  is  a  normalization  constant  and  is  one  of  the  sampling  values  of  a.  For  a 

given  distribution,  the  EUD  and  the  corresponding  figure  of  merit  of  an  IMRT  plan  vary  with  the  sampling 
of  a.  We  thus  rewrite  Eqs.  (3)  and  (4)  as  conditional  probabilities  for  a  sampled  Ok. 

P,{EUD\ak)  =  -^ - - 


^{EUD\ak): 


EUDo 


The  objective  function  for  a  structure  m  in  the  presence  of  uncertainty  in  a  is  expressed  as  the 
summation  of  a  series  of  joint  probabilities 

PSEUD)  =  'Zpseud  I  ak)-PSak) 

*  ’ 

and  the  overall  objective  function  P  of  the  system  is  a  product  of  PJiJEUD)  defined  in  Eq.  (8).  That  is 
F  =  ln(l/P)  =  -\nY[PSEUD)  =  -YS^f^PSEUD  \ 

m  m  k  ^  ^9) 


Figure  1  (Left)  The  target  and  OAR  DVHs  of  four  optimal  plans  when  parameter  a  is  a  fixed  value  (bar  chart  dl)  and  varies 
according  to  three  different  probabilistic  distributions  (bar  chart  d2,  d3  and  d4). 

Figure  2  (Right)  The  EUD  of  the  target  and  objective  Unction  when  parameter  a  is  prescribed  according  to  Fig.  1. 


We  first  investigated  the  behavior  of  the  system  using  a  C-shaped  tumor  case  when  the  parameter  a  of 

the  target  EUD  takes  four  different  distributions,  as  depicted  in 
the  bar  charts  shown  on  the  right  of  Fig.  1,  while  keeping  the 
parameter  a  of  the  OAR  at  a  constant  ao  =  6.0.  In  the  case  shown 
in  Fig.  1  dl,  the  parameter  a  takes  only  a  single  value,  ao  =  -  \Q, 
which  is  a  simple  case  studied  by  Wu  et  al  The  optimal  plans 
for  the  four  distributions  of  parameter  a  differ  significantly,  as 
indicated  by  the  target  and  OAR  DVHs  shown  in  Fig.  1  A  and 
B.  To  estimate  the  degree  of  sensitivity  of  the  solutions  against  a 
variation  in  a,  we  computed  the  target  EUD  and  the  objective 
function,  fr,  as  a  function  of  parameter  a  for  the  four  optimal 
dose  distributions  under  different  types  of  uncertainty 
distributions.  The  results  are  plotted  in  Fig.  2.  The  results 
suggest  that  the  EUD  becomes  much  less  sensitive  to  the 
variation  in  parameter  a  in  the  plans  obtained  with  some  “built- 
in”  distributions  in  parameter  a  (i.e.,  plans  corresponding  to 
Figs.  1  d2  to  d4). 

Four  IMRT  plans  with  different  types  of  pre-assumed 
uncertainties  were  generated  for  a  prostate  tumor  case  (Fig.  3A). 
These  include:  (i)  The  a-parameters  for  both  prostate  target  and 
OARs  are 
restricted  to  A 
single  values 
as  listed  in 

Table  1.  This  plan  serves  as  a  reference  whose  DVHs 
are  shown  in  Figs.  4A-C  as  dotted  curves;  (ii)  Only  the 
a-parameter  of  the  prostate  target  takes  a  range  of 
values,  as  depicted  in  the  right  of  Fig.  4A;  (iii)  Only  the 
a-parameter  of  the  rectum  takes  a  range  of  values,  as 
depicted  in  the  right  of  Fig.  4B;  and  (iv)  The  a- 
parameters  of  both  prostate  target  and  the  rectum  were 
allowed  to  take  a  range  of  values,  as  depicted  in  the 
right  of  Fig.  4C. 

DVHs  for  the  plan  using  parameters  defined  in 
Table  1  are  plotted  with  dotted  curves  and  plans  with 
the  inclusion  of  parameter  uncertainty  are  drawn  with 
solid  curves  (Fig.  4).  When  the  parameter  a  in  target 
EUD  takes  a  Poisson  distribution  as  shown  in  the  bar 
chart  of  Fig.  4A,  prostate  dose  homogeneity  is 
significantly  improved.  The  minimum  dose  increases 
from  55  Gy  to  67  Gy,  and  the  maxim  dose  decreases 
slightly  from  82  Gy  to  80Gy.  However  the  volumes 
receiving  radiation  dose  for  rectum,  bladder  and 


Figure  3.  A  transverse  slice  showing  the 
anatomical  structures  delineated  for  the  prostate 
tumor  (A)  and  the  corresponding  optimized 
dose  distribution  with  the  parameters  listed  in 
Tab.  1  and  the  probabilistic  distribution  shown 
in  Fig.  4B. 
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normal  tissue  all  increase  significantly  though  the 
maximum  dose  remains  similar.  The  improvement  of 
the  target  coverage  and  compromise  of  OAR  sparing 
is  a  natural  outcome  of  the  competitive  requirements 
for  targets  and  OARs  imposed  on  the  system.  The 
corresponding  dose  distribution  with  the  target 


Figure  4.  DVHs  for  a  prostate  cancer  case  using  the 
conventional  optimization  with  fixed  a-value  (dotted  line)  and 
the  newly  proposed  approach  with  the  inclusion  of  model 
parameter  uncertainty  (solid  line).  (A)  Only  the  a-parameter  for 
the  target  is  assigned  with  a  probabilistic  distribution;  (B)  Only 
the  a-parameter  for  the  OAR  is  assigned  with  a  probabilistic 
distribution;  (C)  Uncertainties  in  the  a-parameter  are 
introduced  for  both  the  target  and  OAR. 
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corresponding  dose  distribution  with  the  target  parameter  defined  in  the  bar  chart  A  is  shown  in  Fig.  3B. 
Next  we  considered  the  inclusion  of  parameter  a  uncertainty  in  EUD  calculation  in  one  of  the  critical 
structures-rectum  (Fig.  4  B).  The  irradiated  rectum  volume  for  a  dose  below  60  Gy  is  less  than  that  of  a 
conventional  plan  with  the  parameter  a  fixed  at  24.  DVHs  for  the  bladder,  normal  tissue  and  prostate  do  not 
change  significantly  compared  to  the  plan  without  inclusion  of  parameter  uncertainty.  Lastly,  we 
simultaneously  replaced  target  and  rectum  parameters  with  the  distributions  shown  in  Fig.  4  C.  Similar  to 
that  corresponds  the  prescription  of  Fig.  4  A,  the  prostate  coverage  is  improved.  However,  the  rectum  DVH 
in  this  case  is  not  worsen  greatly  because  parameter  a  of  rectum  EUD  was  allowed  to  take  a  spectrum  of 
values.  For  bladder  and  normal  tissue,  although  their  irradiated  volumes  in  the  low  dose  region  are  higher 
than  those  of  the  conventional  plan,  the  volumes  receiving  high  doses  are  reduced. 

Table  1.  The  conventional  EUD-based  optimization  parameter  for  prostate  cancer. 


PTV 

PTV* 

Bladder 

Rectum 

NT 

a 

-10.0 

10.0 

6.0 

24 

6.0 

EUDo(Gy) 

72 

76 

35 

35 

35 

n 

20 

20 

6 

6 

6 

*  Contains  parameters  for  the  target  treated  as  virtual  normal  tissue  to  limit  dose  inhomogeniety. 


III.  KEY  RESEARCH  ACCOMPLISHMENTS 

•  Developed  a  series  of  mathematical  formulae  to  incorporate  model  parameter  uncertainty  into  IMRT 
optimization. 

•  Verified  the  new  method  in  the  prostate  cancer  Case  and  better  dose  coverage/sparing  can  be  obtained 
with  appropriate  parameters. 


IV.  REPORTABLE  OUTCOMES 

The  following  is  a  list  of  publications  resulted  from  the  grant  support.  Copies  of  the  publication 
materials  are  enclosed  with  this  report. 

Refereed  publication: 

1.  lian  J.,  Cotrutz  C.  and  Xing  L.  Therapeutic  treatment  plan  optimization  with  probability  density-based 
dose  prescription.  Medical  Physics  30:  655-666  (2003). 

2.  lian  J;  and  Xing  L.  Incorporating  Model  Parameter  Uncertainty  into  Inverse  Treatment  Planning, 
Medical  Physics,  (tentatively  accepted),  2004. 

Published  Abstracts: 

1.  Lian  J.,  Spielman  D.,  Cotrutz  C.,  Hunjan  S.,  Adalsteinsson  E.,  King  C.,  Luxton  G.,  Kim  D.,  Daniel  B. 
and  Xing  L.  Including  metabolic  uncertainty  into  proton  MR  spectroscopic  imaging  (MRSI)-guided  inverse 
treatment  planning.  In:  the  45th  American  Association  of  Physicists  in  Medicine  (AAPM)  Meeting,  San 
Diego,  2003 

2.  Lian  J.  and  Xing  L.  Biological  Model  Based  IMRT  Optimization  with  Inclusion  of  Parameter 
Uncertainty.  In;  the  14th  International  Conference  on  the  Use  of  Computers  in  Radiation  Therapy  (ICCR), 
Seoul,  Korea,  2004 

New  Employment: 

With  the  help  of  this  grant  and  associated  publications,  PI  has  obtained  an  assistant  professor  position 
in  the  Department  of  Radiation  Oncology  at  the  University  of  North  Carolina  at  Chapel  Hill. 
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V.  CONCLUSIONS 


A  technique  for  incorporating  biological  model  parameter  uncertainties  into  inverse  treatment  planning 
has  been  developed.  The  formalism  is  quite  general  and  does  not  prerequisite  the  specific  form  of 
uncertainty  distributions  of  the  involved  model  parameters.  By  including  model  parameter  uncertainties,  the 
final  solution  becomes  more  robust  and  the  treatment  outcome  will  be  less  likely  influenced  by  inter-patient 
variation  of  biological  characteristics.  With  the  increasing  interest  in  radiation  therapy  community  to  use 
biologically  based  models  for  treatment  planning,  this  work  provides  an  effective  way  to  better  account  for 
the  known  uncertainties  in  the  model  parameters  and  allows  us  to  maximally  utilize  the  available 
radiobiology  knowledge  to  facilitate  patient  care. 
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Therapeutic  treatment  plan  optimization  with  probability  density-based 
dose  prescription 

Jun  Lian,  Cristian  Cotrutz,  and  Lei  Xing®^ 

Department  of  Radiation  Oncology,  Stanford  University  School  of  Medicine,  300  Pasteur  Drive,  Stanford, 

California  94305-5304 

(Received  18  September  2002;  accepted  for  publication  29  January  2003;  published  26  March  2003) 

The  dose  optimization  in  inverse  planning  is  realized  under  the  guidance  of  an  objective  function. 

The  prescription  doses  in  a  conventional  approach  are  usually  rigid  values,  defining  in  most  in¬ 
stances  an  ill-conditioned  optimization  problem.  In  this  work,  we  propose  a  more  general  dose 
optimization  scheme  based  on  a  statistical  formalism  [Xing  et  al,  Med.  Phys.  21,  2348-2358 
(1999)].  Instead  of  a  rigid  dose,  the  prescription  to  a  structure  is  specified  by  a  preference  function, 
which  describes  the  user’s  preference  over  other  doses  in  case  the  most  desired  dose  is  not  attain¬ 
able.  The  variation  range  of  the  prescription  dose  and  the  shape  of  the  preference  function  are 
predesigned  by  the  user  based  on  prior  clinical  experience.  Consequently,  during  the  iterative 
optimization  process,  the  prescription  dose  is  allowed  to  deviate,  with  a  certain  preference  level, 
from  the  most  desired  dose.  By  not  restricting  the  prescription  dose  to  a  fixed  value,  the  optimiza¬ 
tion  problem  becomes  less  ill-defined.  The  conventional  inverse  planning  algorithm  represents  a 
special  case  of  the  new  formalism.  An  iterative  dose  optimization  algorithm  is  used  to  optimize  the 
system.  The  performance  of  the  proposed  technique  is  systematically  studied  using  a  hypothetical 
C-shaped  tumor  with  an  abutting  circular  critical  structure  and  a  prostate  case.  It  is  shown  that  the 
final  dose  distribution  can  be  manipulated  flexibly  by  tuning  the  shape  of  the  preference  function 
and  that  using  a  preference  function  can  lead  to  optimized  dose  distributions  in  accordance  with  the 
planner’s  specification.  The  proposed  framework  offers  an  effective  mechanism  to  formalize  the 
planner’s  priorities  over  different  possible  clinical  scenarios  and  incorporate  them  into  dose  opti¬ 
mization.  The  enhanced  control  over  the  final  plan  may  greatly  facilitate  the  IMRT  treatment 
planning  process.  ©  2003  American  Association  of  Physicists  in  Medicine. 

[DOI:  10.1118/1.1561622] 
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I.  INTRODUCTION 

Inverse  planning  is  used  in  intensity  modulated  radiation 
therapy  (IMRT)  for  deriving  the  optimal  beam  intensity  pro¬ 
files  that  produce  the  best  possible  dose  distribution  for  a 
given  patient. The  dose  optimization  process  is  usually 
performed  under  the  guidance  of  an  objective  function, 
which  measures  the  “distance”  between  the  physical  and  the 
prescribed  dose  distributions.^’^^"^®  One  of  the  common  ob¬ 
jective  functions  for  inverse  planning  is  the  quadratic  objec¬ 
tive  function,^’^^’^^  with  importance  factors  assigned  to  the 
involved  structures  to  prioritize  their  relative  importance  dur¬ 
ing  the  optimization  process.^^"^^  The  objective  function  is 
defined  as  a  global  quantity  based  on  general  physical  con¬ 
siderations.  When  the  desired  dose  distribution  is  not  attain¬ 
able  during  optimization,  a  compromise  solution  is  found 
using  the  algorithm’s  ranking.  The  compromise  dose  distri¬ 
bution,  however,  is  often  not  what  the  planner  wants  and 
multiple  trial  and  errors  are  needed  to  obtain  a  clinically 
acceptable  IMRT  plan. 

A  main  problem  of  the  existing  IMRT  planning  algo¬ 
rithms  is  the  lack  of  an  effective  mechanism  for  incorporat¬ 
ing  prior  knowledge  into  inverse  planning.^*  In  the  past, 
there  have  been  many  attempts  to  introduce  soft/hard  con¬ 
straints  to  steer  the  dose  optimization  process  toward  the 


clinically  desired  solutions. However,  the  constraints 
are  introduced  in  an  ad  hoc  fashion  and  do  not  fully  utilize 
the  partial  information  available  from  years  of  clinical  inves¬ 
tigations  because  of  their  phenomenological  nature.  On  a 
more  fundamental  level,  the  constraints  are  imposed  a  pos¬ 
teriori  and  controls  the  optimization  passively.  Our  purpose 
in  this  paper  is  to  develop  a  statistical  analysis-based  inverse 
planning  formalism  to  more  effectively  utilize  the  prior 
knowledge.  Instead  of  specifying  a  rigid  prescription  dose, 
the  formalism  allows  us  to  use  a  dose  distribution  as  the 
input  prescription  to  the  system,  providing  a  natural  way  for 
us  to  take  advantage  of  the  existing  information  of  the  sys¬ 
tem  variables  and  promising  to  make  the  optimization  out¬ 
come  more  predictable  and  controllable. 

In  the  next  section  we  present  the  details  of  the  new  dose 
optimization  algorithm  after  a  brief  introduction  of  the  con¬ 
cept  of  preference  function.  The  formalism  is  then  applied  to 
a  synthetic  phantom  case  with  C-shaped  tumor  target  and  a 
prostate  case.  Our  results  indicate  that  the  statistical  analysis- 
based  formalism  provides  a  general  framework  for  inverse 
planning  and  is  capable  of  producing  conformal  IMRT  dose 
distribution.  Coupled  with  the  capability  of  the  preference 
function  in  customizing/formalizing  our  prior  clinical  knowl¬ 
edge,  it  is  expected  that  the  proposed  technique  will  have  a 
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broad  implication  and  potential  to  greatly  facilitate  an  IMRT 
planning  process. 

11.  MATERIAL  AND  METHODS 

A.  Theoretical  background 

In  a  vectorial  form,  the  dose  to  the  points  in  the  treatment 
region  depend  upon  the  beamlet  weights  w  as 

Dc  =  d-w,  (1) 

where  d  represents  the  dose  deposition  matrix,  expressing 
the  dose  deposited  to  any  point  in  the  patient  when  irradiated 
with  a  unit  weight  beamlet  vector.  The  inverse  problem  as 
posed  for  IMRT  is  to  find  a  set  of  beamlet  weights  that  pro¬ 
duce  the  optimal  dose  distribution  by  minimizing  a  therapeu¬ 
tic  objective  function.  The  most  used  objective  function  has  a 
quadratic  form  and  reads^^  as 

1 

ra[D,(n)-Dpin)]\  (2) 

where  N  is  the  total  number  of  voxels,  is  the  importance 
factor  that  controls  the  relative  importance  of  a  structure  cr, 
and  Dp  and  are  prescribed  and  calculated  doses,  respec¬ 
tively. 

In  inverse  planning  algorithm  based  on  the  quadratic  ob¬ 
jective  function  [Eq.  (2)],  the  dose  prescription  to  the  target 
or  sensitive  structure  takes  a  rigid  value.  The  minimization  of 
the  objective  function  is  realized  by  various  algorithms  like 
simulated  annealing,  gradient  methods,  Indepen¬ 

dent  of  the  used  dose  optimization  algorithms,  we  will  call 
these  methods  throughout  the  text  conventional  IMRT  opti¬ 
mization  procedures.  The  problem  is  usually  ill-posed  and 
may  lead  to  negative  fluence  unless  hard  constraints  are 
introduced.^  Practically,  it  is  not  uncommon  that  the  plans 
computed  by  what  are  called  optimization  systems  are  not 
consistent  with  the  expectation  of  the  planner  and  that  sev¬ 
eral  trial-and-error  adjustments  of  the  system  parameters 
might  be  required  to  achieve  a  clinically  acceptable  plan. 
Given  a  patient,  the  obtained  plan  can  vary  widely  from  one 
planer  to  the  next,  even  within  a  department.  In  the  following 
we  describe  a  more  adaptable  and  “intelligent”  statistical 
inverse  planning  formalism  based  on  the  concept  of  a  pref¬ 
erence  function  to  better  deal  with  the  dilemma. 

B.  Preference  function 

In  a  recent  paper,  Xing  et  al?^  introduced  the  concept  of 
preference  function  to  weaken  the  rigid  dose  prescription 
commonly  seen  in  the  existing  inverse  planning  algorithms. 
Its  role  is  to  allow  a  dose  distribution  to  be  considered  in¬ 
stead  of  just  a  single  value,  and  to  quantify  the  degree  of  our 
willingness  to  accept  a  prescription  dose  Dp  in  that  range. 
The  preference  function  can  be  constructed  heuristically 
from  clinical  considerations.^^  The  defined  preference  func¬ 
tion  states  that  the  most  favorable  prescription  dose  for  a 
voxel  n\sDp{n)  and  that  a  different  prescription  dose  is  also 
acceptable,  but  with  a  smaller  preference  level.  For  illustra¬ 
tion,  in  Fig.  1  we  show  a  sketch  of  the  preference  functions 
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Fig.  1.  A  sketch  of  preference  functions  for  a  target  and  a  sensitive  structure. 


for  a  target  and  sensitive  structure.  The  most  desirable  dose 
for  a  sensitive  structure  should  generally  be  set  to  zero.  The 
conventional  prescription  scheme  represents  a  special  case  of 
the  general  approach  proposed  here  with  the  step  function 
form  of  the  preference  function.  That  is. 


PniDp)^ 


1, 

0, 


if  D,=dI, 

if  Dpi^Dl. 


(3) 


To  give  another  example,  we  write  down  the  Gaussian 
preference  function  for  a  voxel  n 


Pn{Dp)  =  P^n  exp{-  ri£>p(n)-Z)°(n)]2},  (4) 


where  Fq/i  is  a  normalization  constant  and  represents  the 
Gaussian  parameter.  For  a  system  comprising  N  voxels,  the 
total  preference  is  given  by  a  product  of  the  preference  func¬ 
tions  of  all  voxels: 


f’  =  n  PniD,) 

n 


=  n  f’o«exp{-r„[Dp(n)-Z)°(n)]2}.  (5) 

When  a  maximum  likelihood  estimator  is  used,  it  has 
been  demonstrated  that  the  maximization  of  the  logarithmic 
function  of  P  or  minimization  of  ln(l/P),  is  equivalent  to  the 
minimization  of  the  conventional  quadratic  objective 
function.^'’^'*  In  this  case,  the  Gaussian  parameter  in  Eq. 
(5),  which  commands  the  “spread”  of  the  Gaussian  around 
£)p,  is  equivalent  to  the  importance  factor  that  controls  the 
relative  importance  of  the  structure  and  parametrizes  the 
clinical  trade-off  strategy. 

C.  Probability  density-based  dose  prescription  and 
inverse  planning 

The  objective  function  defined  in  Eq.  (2)  uses  a  rigid 
dose.  Dp .  Since  in  most  instances  an  ideal  dose  prescription 
is  not  physically  attainable,  we  resort  to  an  expansion  of  the 
prescription  dose,  over  a  certain  interval.  That  is,  we  allow 
the  prescription  dose  to  take  a  “probabilistic”  distribution 
around  the  most  desired  dose  as  specified  by  the  preference 
function.  For  computational  purpose,  we  divide  the  permis- 
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sible  prescription  dose  into  a  number  of  discretized  values, 
where  i  is  the  index  of  a  possible  prescription  dose  and 
f  =  0  represents  the  most  desirable  dose.  The  preference  dis¬ 
tribution  prescription  is  usually  normalized  to  unity. 

In  order  to  utilize  the  probability  information  character¬ 
ized  by  the  preference  function,  we  formulate  the  conven¬ 
tional  dose  optimization  into  a  statistical  analysis  problem. 
To  proceed,  let  us  take  the  quadratic  objective  function  as  an 
example.  We  rewrite  the  traditional  quadratic  objective  func¬ 
tion  (2)  into 

/(^c)=/on  exp{-r„[Z)^(n)-I>  (rt)]^}.  (6) 

n  ^ 

where  /o  is  a  normalization  constant.  For  a  given  prescribed 
dose  distribution,  Eq.  (6)  measures  the  goodness  of  a  calcu¬ 
lated  dose  distribution  using  an  exponential  scale,  as  com¬ 
pared  with  Eq.  (2).  Equation  (6)  can  be  interpreted  as  a  con¬ 
ditional  probability  and  formally  rewritten  as 

fiDc\Dp)=fQj[  exp{-r<,[£)^(n)-Z)p(n)]2}.  (7) 

n 

When  the  prescription  dose  is  no  longer  a  rigid  dose,  it  is 
conceivable  that  there  are  a  number  of  optimum  solutions, 
each  corresponding  to  a  sample  of  prescription  doses.  Math¬ 
ematically,  we  now  have  two  “probability”  distribution  func¬ 
tions.  One  is  the  preference  function  that  characterizes  our  a 
priori  preference  over  different  prescription  doses  P(Dp), 
and  the  other  is  Eq.  (6)  that  ranks  a  calculated  dose  for  a 
given  prescribed  dose.  Dp .  Our  task  is  to  find  the  solution 
that  is  statistically  optimal  with  consideration  of  the  variable 
prescription.  For  this  purpose,  we  introduce  the  “joint  prob¬ 
ability”  of  the  two  “probability”  distributions  defined  by 
Eqs.  (5)  and  (7).  The  function  at  a  voxel  n  can  be  written  as 

^(^c)  =  2  /„(Z),|D;)P„(Dp.  (8) 

The  total  preference  function  of  the  system  is  given  by 

/’=n  Pn(D,).  (9) 

n 

D.  Optimization  strategy 

Having  the  rigid  prescription  Dp  in  (2)  replaced  by  a 
range  of  prescribed  doses,  {7)^},  the  total  preference  func¬ 
tion  is  now  given  by  Eqs.  (8)  and  (9).  For  convenience,  we 
define  objective  function  F=ln(l/F)  and  derive  the  optimal 
solution  by  minimizing  the  F,  which  is  equivalent  to  maxi¬ 
mize  the  preference  function  (9).  The  objective  function  now 
reads  as 

F=ln(l/P) 

=  -lnn^’„=-S  lnE/„(0,|D;)  P„(Z);).  (10) 

n  n  i  y  y 


Fro.  2.  A  flow  chart  of  the  optimization  process  with  the  inclusion  of  pre¬ 
designed  preference  function  information. 


Note  that  the  conventional  quadratic  objective  function  is  a 
special  case  of  the  above  general  objective  function  when  the 
prescription  takes  a  rigid  value  for  each  structure,  as  de¬ 
scribed  by  Eq.  (3). 

The  optimization  process  is  schematically  shown  in  a  flow 
chart  (Fig.  2).  The  beam  profile  is  determined  by  minimizing 
the  above  objective  function  using  a  conjugate  gradient  op¬ 
timization  algorithm.  The  details  of  the  algorithm  have  been 
discussed  in  a  previous  paper.^^  Briefly,  the  calculation  con¬ 
sists  of  three  major  steps:  (i)  assume  an  initial  intensity  pro¬ 
file  for  each  incident  beam;  (ii)  compute  the  “joint  probabil¬ 
ity”  given  by  Eqs.  (8)  and  (9).  For  this  purpose,  we  need  to 
sample  all  combinations  of  the  prescription  doses  of  different 
structures  and  compute  the  function  given  in  Eqs.  (5)  and  (7) 
for  each  of  these  combinations;  and  (iii)  optimization  of  the 
multidimensional  “joint  probability”  function.  The  second 
step  is  fairly  computationally  intensive  because  we  must 
compute  the  two  functions  for  every  sampling  of  the  pre¬ 
scription  doses.  In  our  calculation,  we  typically  assign  four 
to  seven  discrete  possible  prescription  doses  for  each  struc¬ 
ture.  A  finer  discretization  of  the  prescription  dose  did  not 
seem  lead  to  further  improvement  but  would  greatly  increase 
the  computation  time.  All  calculations  presented  here  are 
performed  on  a  Personal  Computer  (PC)  with  an  Intel  Pen¬ 
tium®  in  1  GHz  CPU  (Intel  Corporation,  Sunnyvale,  CA). 
The  computation  time  needed  to  obtain  an  optimal  solution 
for  a  given  set  of  system  parameters  (including  beam  con¬ 
figuration,  preference  function,  importance  factors)  is  typi¬ 
cally  less  than  ten  minutes. 

III.  RESULTS  AND  DISCUSSION 

A.  A  synthetic  phantom  case  with  a  C-shaped  tumor 

To  systematically  study  the  performance  of  the  statistical 
analysis-based  inverse  planning  algorithm,  we  applied  the 
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A 


B 


Fig.  3.  (a)  A  sketch  of  a  phantom  case  with  a  C-shaped 
tumor.  The  dose  prescription  is  set  100  (arbitrary  units)  to 
the  FTV  and  0  to  the  circular  OAR  and  normal  tissue,  (b) 
Dose  distribution  obtained  using  the  “probabilistic”  pre¬ 
scription  shown  in  Fig.  4(a). 


technique  to  a  C-shaped  tumor  case  [Fig.  3(a)]  with  a  variety 
of  preference  functions  and  compared  the  results  with  that 
obtained  using  the  conventional  approach  with  a  fixed  dose 
prescription.  Nine  equally  spaced  6  MV  beams  beginning  at 
Qo  ^ggJ  study.  The  prescription  doses  to 

the  PTV  and  OAR  in  the  conventional  IMRT  plans  were  100 
and  0  (the  dose  is  in  an  arbitrary  unit),  respectively. 

We  first  assigned  three  sets  of  symmetrical  Gaussian  dis¬ 
tributions  to  the  target  while  keeping  the  prescription  to  the 
sensitive  structure  at  zero  (Fig.  4).  The  Gaussian  preference 
functions  were  represented  by  three  sets  of  preference  levels 
at  seven  discrete  values  (80,  87,  94,  100,  106,  113,  and  120). 
The  center  of  the  Gaussian  functions  was  set  at  100.  The 
preference  levels  for  the  seven  doses  are  shown  in  Fig.  4  for 
each  of  the  three  situations  studied  here.  The  transverse  dose 
distribution  obtained  using  the  statistical  inverse  planning 
formalism  for  the  case  shown  in  Fig.  4(a)  is  plotted  in  Fig. 
3(b).  As  expected,  target  inhomogeneity  increases  as  we 


loosen  the  constraint  of  the  rigid  dose  prescription.  This  can 
be  better  demonstrated  by  using  the  differential  DVH  for 
each  situation.  As  seen  from  the  differential  DVH  plots  (the 
right  column  of  Fig.  4),  the  width  of  the  differential  function 
gradually  increases,  from  26.72,  28.59-30.39,  as  we  gradu¬ 
ally  increase  the  acceptance  levels  for  the  doses  different 
from  the  most  desirable  dose  (100).  This  series  of  calcula¬ 
tions  provides  us  with  preliminary  evidence  that  the  final 
dose  distribution  can  be  steered  by  varying  the  preference 
function. 

Next,  we  constructed  six  sets  of  asymmetric  preference 
functions  for  the  target  (Fig.  5  and  Fig.  6).  When  higher 
preference  levels  were  assigned  to  the  doses  higher  than  100, 
we  found  that  the  target  DVH  is  shifted  to  the  high  dose 
region.  Interestingly,  even  when  an  extremely  low  preference 
(for  instance,  1%)  was  assigned  to  the  doses  less  than  100 
[Fig.  5(b)],  a  noticeable  underdosing  relative  to  the  conven¬ 
tional  result  was  resulted.  A  similar  phenomenon  can  also  be 
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Pose  {^bitrary  unit) 


Fig.  4.  DVHs  of  the  PTV,  OAR,  and  normal  tissue  (NT)  obtained  using  the  conventional  rigid  dose  prescription  (dotted  line)  and  the  “probabilistic” 
prescription  (solid  line).  The  Gaussian  preference  functions  with  different  variances  are  shown  in  the  middle  panel.  The  right  panel  shows  the  differential 
DVHs  for  the  target. 
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Dose  (arbitrary  unit) 
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Dose  (arbitrary  unit) 


Dpse(arbltraivvinit) 
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Rg.  5.  DVHs  of  the  PTV,  OAR,  and  normal  tissue  (NT)  obtained  with  the  conventional  rigid  dose  prescription  (dotted  line)  and  with  the  “probabilistic” 
prescription  (solid  line).  The  bar  charts  on  right  show  the  asymmetrical  preference  functions. 
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Dose  (arbitri^  iihit) 


Dose  (arbitrary  unit) 


Rg.  6.  DVHs  of  the  PTV,  OAR,  and  normal  tissue  (NT)  obtained  using  the  conventional  rigid  dose  prescription  (dotted  line)  and  the  new  statistical  inverse 
planning  method  for  a  variety  of  preference  functions  shown  on  the  right  panel  (solid  line). 


seen  from  the  result  shown  in  Fig.  5(c),  where  only  0.5%, 
0.7%,  and  1%  of  preference  levels  were  assigned  to  the  dose 
values  of  80,  87,  and  94.  This  observation  seems  to  indicate 
that  the  influence  of  the  assigned  preference  level  at  a  low 


dose  plays  an  important  role.  In  Fig.  6,  we  set  the  preference 
levels  for  the  doses  less  than  100  to  be  0  and  only  assign 
nonzero  preference  levels  for  the  doses  higher  than  100.  It  is 
seen  that  in  all  these  situations  the  minimum  target  dose  is 
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Dose  tiii^bltrgry  unit) 


Fig.  7.  DVHs  of  the  OAR  and  PTV  when  the  prescription  dose  to  the  OAR  is  modeled  by  (a)  a  uniform  distribution,  (b)  a  bell-shaped  function,  (c)  an 
exponential  decay  function,  and  (d)  a  rigid  value. 

higher  than  that  of  the  conventional  plan.  As  a  result  of  our  OAR  dose  to  take  seven  values:  0,  5,  10,  15,  20,  25,  and  30 

preference  over  higher  doses,  the  fractional  volume  at  any  with  the  acceptance  levels  sampled  from  three  different  types 

dose  less  than  100  is  improved  in  comparison  to  that  of  the  of  prescription  distribution:  uniform  [Fig.  7(a)],  bell-shaped 

conventional  IMRT  plan.  In  Fig.  6(c),  we  further  exemplify  [Fig.  7(b)],  and  exponential  [Fig.  7(c)]  functions.  Figure  7(d) 

the  statistical  analysis  based  inverse  planning  method  by  represents  the  conventional  case  with  zero  prescription  to  the 

simplifying  our  preference  to  two  doses  (100  and  120),  each  OAR.  The  corresponding  OAR  and  PTV  DVHs  are  plotted 

with  50%  preference  levels.  In  this  situation,  in  addition  to  in  the  left  panel  of  Fig.  7.  When  the  preference  was  uni- 

that  the  doses  in  the  target  are  shifted  toward  higher  values,  formly  sampled  in  the  dose  interval  from  0  to  30,  the  result- 

the  target  DVH  exhibits  a  stepwise  behavior:  a  plateau  ap-  ant  dose  to  the  OAR  was  found  to  be  the  highest,  as  indi¬ 

pears  at  around  110,  which  is  in  the  middle  of  the  two  pre-  cated  by  curve  A  in  Fig.  7.  The  best  target  dose  coverage  was 
scribed  doses.  achieved  in  this  situation.  If  the  preference  to  a  high  dose 

It  is  interesting  to  point  out  that  the  OAR  sparing  is  im-  was  reduced,  the  DVH  was  gradually  shifted  to  the  low  dose 

proved  as  compared  with  the  conventional  IMRT  plan  in  direction  (curves  B  and  C).  It  is  not  surprising  that  the  best 

most  cases  studied  in  Figs.  5  and  6,  even  when  the  target  OAR  sparing  was  achieved  in  the  conventional  case  where  a 

dose  is  escalated.  That  is,  the  DVH  of  the  OAR  is  not  always  zero  dose  was  prescribed  to  the  OAR.  The  target  dose  homo- 

shifted  toward  higher  doses,  as  would  occur  if  a  higher  dose  geneity  was  slightly  improved  in  all  cases  when  a  probabi- 

is  prescribed  in  a  conventional  inverse  planning  system.  In-  listic  prescription  was  given  to  the  OAR.  Similar  to  that  de- 

stead,  the  dose  to  the  OAR  remained  unchanged  or  even  scribed  in  the  last  paragraph,  the  results  clearly  demonstrate 

lowered  in  some  cases.  A  reasonable  explanation  for  the  ob-  that  the  “probabilistic”  prescription  allows  us  to  control  the 

served  phenonienon  is  that,  when  a  rigid  dose  prescription  is  OAR  dose  distribution  and  indicate  the  usefulness  of  the 

replaced  by  a  range  of  doses,  the  system  is  given  more  free-  statistic  analysis  approach, 

dom  for  self-adjustment.  As  a  benefit,  a  solution  with  a 
higher  integral  target  dose  and  reduced  OAR  dose  can  be 

obtained  from  the  expanded  solution  space.  ’  ®  prostate  case 

We  have  also  studied  the  behavior  of  the  system  when  a  The  new  inverse  planning  algorithm  was  also  applied  to 
range  of  doses  is  prescribed  to  the  OAR.  In  this  investiga-  study  a  six  filed  IMRT  prostate  treatment  [Fig.  8(a)].  Four 

tion,  we  kept  the  target  prescription  to  100  and  allowed  the  plans  with  different  types  of  preference  functions  were  gen- 
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Fig.  8.  A  transverse  slice  showing  the  anatomical  struc¬ 
tures  delineated  for  the  prostate  tumor  (a)  and  (b)  the 
dose  distribution  obtained  using  the  “probabilistic” 
prescription  shown  in  Fig.  9(a). 


crated.  In  addition,  a  plan  with  rigid  prescription  (74  Gy  on 
the  target,  60  Gy  on  the  bladder,  and  40  Gy  on  the  rectum)  is 
also  generated.  The  DVHs  for  this  plan  is  plotted  as  dotted 
lines  in  Fig.  9  and  is  used  as  a  reference  for  comparison.  In 
all  treatment  plans,  six  beams  were  placed  at  the  following 
angular  positions:  0°,  55°,  135°,  180°,  225°,  and  305°.  The 
size  of  the  pencil  beam  defined  at  the  isocenter  was  0.5  cm. 
The  DVH  and  preference  functions  for  four  different 
plans  are  schematically  shown  in  Fig.  9.  In  the  study  shown 
in  Figs.  9(a)“9(b),  we  kept  the  preference  function  of  the 
sensitive  structures  unchanged  and  only  varied  the  form  of 
the  preference  function  of  the  target.  In  Fig.  9(a),  we  as- 
sumed  that  target  could  take  seven  discrete  values  (74,  76, 
78,  80,  82,  84,  and  86  Gy)  sampled  from  an  exponential 
distribution.  Compared  with  the  dotted  lines,  the  target  DVH 
was  shifted  toward  the  high  dose  direction.  The  dose  distri¬ 
bution  corresponding  to  the  preference  function  is  shown  in 
Fig.  8(b),  The  target  DVH  was  shifted  even  further  toward 


the  high  dose  region  [Fig.  9(b)]  when  a  bell-shaped  prefer¬ 
ence  function  was  used  with  more  emphasis  on  the  target 
receiving  doses  at  74,  76,  and  78  Gy.  In  both  cases,  doses  to 
the  rectum  and  bladder  did  not  change  significantly. 

In  Fig.  9(c)  we  show  the  DVHs  when  the  preference  func¬ 
tion  to  the  rectum  deviates  from  the  uniform  distribution.  As 
a  result,  the  rectum  dose  was  significantly  lowered  in  all  dose 
levels  and  the  maximum  dose  was  reduced  from  66  to  57  Gy. 
Because  of  the  proximity  of  the  rectum  to  the  prostate  target, 
the  maximum  rectum  dose  was  not  restricted  to  30  Gy,  as 
specified  in  the  preference  function.  We  emphasize  that  the 
improvement  in  rectum  and  bladder  sparing  was  achieved  at 
cost  of  higher  dose  inhomogeneity  in  the  prostate  target.  This 
reminds  us  that,  in  dose  optimization,  there  is  a  dosimetric 
compromise.  That  is,  the  improvement  in  the  dose  to  a  struc¬ 
ture  is  often  accompanied  by  dosimetrically  adverse  effect(s) 
at  other  points  in  the  same  or  different  structures.  The  impor¬ 
tant  point  that  one  should  note  is  that  from  the  clinical  point 
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Fig.  9.  A  comparison  of  the  DVHs  obtained  using  the  “probabilistic”  prescription  (solid  line)  and  conventional  rigid  dose  prescription  (dotted  line).  The 
prostate  and  rectum  prescriptions  are  represented  on  the  right  bar  charts. 
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of  view,  some  dose  distributions  are  more  acceptable  than 
others  and  our  goal  is  to  find  the  solution  that  improves  the 
plan  to  the  largest  possible  extent,  but  with  a  clinically  insig¬ 
nificant  or  acceptable  sacrifice  in  OAR  sparing.  In  order  to 
achieve  this,  it  is  necessary  to  have  a  reasonable  amount  of 
controllability  over  the  final  dose  distribution.  In  this  sense 
we  believe  that  the  proposed  formalism  is  valuable. 

In  addition,  we  varied  the  preference  functions  for  both 
target  and  rectum  [Fig.  9(d)].  The  preference  function  for  the 
rectum  was  the  same  as  that  in  Fig.  9(c).  Compared  to  the 
results  shown  in  Fig.  9(c),  we  found  that  the  target  dose 
inhomogeneity  was  slightly  improved. 


IV.  CONCLUSIONS 

The  formalism  we  derived  here  provides  a  general  starting 
point  for  the  study  of  a  system  with  a  probability  density- 
based  dose  prescription.  The  inclusion  of  the  partial  informa¬ 
tion  into  the  plan  selection  process  represents  a  significant 
change  from  the  conventional  approaches.  The  proposed 
technique  can  be  categorized  into  the  general  Bayesian 
decision-making  theory,^  ^  which  is  a  useful  tool  to  deal  with 
a  system  with  “statistical*’  inference.  In  image  analysis  and 
many  other  fields  of  science  and  engineering,  it  has  proven 
extremely  useful  to  include  the  prior  knowledge  of  the  sys¬ 
tem  variables  into  the  estimation  process.^ ^  The  preference 
function  proposed  for  radiotherapy  optimization  here  serves 
as  a  priori  probability  density  function  in  standard  Bayesian 
statistics.  The  role  of  the  preference  function  is  to  indicate 
our  “bias”  on  the  values  of  the  system  variables.  By  utilizing 
the  partial  information  of  the  system  variables,  one  can  more 
effectively  search  the  solution  space  and  eliminate  some  un¬ 
necessary  uncertainties  in  the  optimization  process. 

In  conclusion,  we  have  developed  a  statistical  analysis- 
based  inverse  planning  algorithm  to  include  preference  and 
expert  knowledge  into  the  dose  optimization  process.  Instead 
of  a  rigid  dose  prescription,  the  new  approach  allows  us  to 
prescribe  a  range  of  doses  with  predesigned  preference  lev¬ 
els.  The  technique  represents  a  novel  application  of  the  gen¬ 
eral  Bayesian  decision-making  theory^ ^  for  dealing  with  sta¬ 
tistical  inference  and  is  valuable  for  deriving  a  statistically 
optimal  solution  in  the  presence  of  uncertainties  in  system 
parameters.  The  method  was  demonstrated  for  a  system  with 
modulating  prescriptions  but  can  be  easily  extended  to  solve 
many  other  related  problems  (e.g.,  in  biologically  based  dose 
optimization,  one  can  incorporate  the  uncertainties  of  various 
radiobiology  parameters  into  the  inverse  planning  process 
using  the  frameset  developed  in  this  work^^).  The  ill  condi¬ 
tioning  of  the  problem  was  improved  because  of  the  use  of  a 
less  restrictive  prescription  and,  as  a  result,  new  solutions 
that  are  otherwise  inaccessible  can  be  obtained  naturally.  It  is 
demonstrated  that  the  obtained  solutions  using  the  new  ap¬ 
proach  strongly  correlate  with  the  preference  function,  sug¬ 
gesting  that  the  planning  process  is  controllable  and  predict¬ 
able  by  the  proposed  method. 
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ABSTRACT 


Radiobiological  treatment  planning  depends  not  only  on  the  accuracy  of  the  models 
describing  the  dose-response  relation  of  different  tumors  and  normal  tissues  but  also  on  the 
accuracy  of  tissue  specific  radiobiological  parameters  in  these  models.  Whereas  the  general 
formalism  remains  the  same,  different  sets  of  model  parameters  lead  to  different  solutions 
and  thus  critically  determine  the  final  plan.  Here  we  describe  an  inverse  planning  formalism 
with  inclusion  of  model  parameter  uncertainties.  This  is  made  possible  by  using  a  statistical 
analysis-based  fi-ameset  developed  by  our  group.  In  this  formalism,  the  uncertainties  of 
model  parameters,  such  as  the  parameter  a  that  describes  tissue-specific  effect  in  EUD 
model,  are  expressed  by  probability  density  function  and  are  included  in  the  dose 
optimization  process.  We  found  that  the  final  solution  strongly  depends  on  distribution 
functions  of  the  model  parameters.  Considering  that  currently  available  models  for 
computing  biological  effects  of  radiation  are  simplistic,  and  the  clinical  data  used  to  derive 
the  models  are  sparse  and  of  questionable  quality,  the  proposed  technique  provides  us  with 
an  effective  tool  to  minimize  the  effect  caused  by  the  uncertainties  in  a  statistical  sense.  With 
the  incorporation  of  the  uncertainties,  the  technique  has  potential  for  us  to  maximally  utilize 
the  available  radiobiology  knowledge  for  better  IMRT  treatment. 

Key  words;  inverse  planning,  dose  optimization,  biological  models,  IMRT 


2 


INTRODUCTION 


Most  IMRT  optimization  systems  at  present  use  dose  and/or  dose  volume-based 
objective  functions  which  guide  the  IMRT  planning  by  imposing  a  penalty  according 
to  the  difference  between  the  computed  and  prescribed  doses.  A  well-known  drawback  of 
the  dose-based  inverse  planning  is  that  the  nonlinear  dose  response  of  tumor  or  normal 
structures  is  not  fully  considered.  A  number  of  mathematical  models  have  been 
developed  over  the  years  to  better  describe  the  biological  effect  of  radiation,  which 
include  tumor  control  probability  (TCP)  normal  tissue  complication  probability 
(NTCP)  equivalent  uniform  dose  (EUD)  ^  and  the  probability  of  uncomplicated  tumor 
control  (P+)  In  parallel  to  these  modeling  efforts,  considerable  works  have  also 

been  done  to  use  these  biological  models  to  construct  more  meaningful  objective 
functions  for  therapeutic  dose  optimization  . 

Generally  speaking,  radiobiological  formalism  involves  the  use  of  model 
parameters  that  are  of  considerable  uncertainty  For  instance,  the  radiosensitivity  a 
of  Webb’s  TCP  model  varies  from  0.157  Gy"'  to  0.090  Gy"'  when  model  parameters  were 
fit  to  103  patients’  data  Biological  ‘margins’  have  been  used  to  account  for  the 
variability  in  radiation  sensitivity.  Similar  to  the  use  of  a  safety  margin  to  account  for  the 
potential  uncertainties,  in  targeting  a  tumor,  this  method  assigns  more  conservative 
radiosensitivity  values  to  the  tumor  or  sensitive  structures  to  deal  with  the  potential 
uncertainty  of  the  parameter^^  Kaver  et  al  proposed  a  stochastic  optimization  to  account 
for  clinical  uncertainties,  including  the  varying  radiosensitivity  The  objective 

function  was  constructed  based  on  a  linear  quadratic  Poisson  model  which  approximates 
the  probability  of  curing  the  patient  or  inflicting  injxuy.  Two  parameters  in  the  model 
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could  be  calculated  if  the  standard  deviation  of  dose  per  faction  was  known.  The 
optimization  was  thus  executed  corresponding  to  different  standard  deviations. 

We  have  recently  presented  a  general  statistical  analysis-based  inverse  planning 
framework  and  applied  it  to  investigate  the  influence  of  model  parameter 

uncertainties  in  biologically  based  dose  optimization  The  purpose  of  this  paper  is  to 
provide  a  detailed  description  of  the  .  technique  and  addresses  several  important  issues 
related  to  the  dose  optimization  in  the  presence  of  model  parameter  uncertainties.  In  our 
approach,  the  uncertainty  of  a  model  parameter  is  quantified  by  a  probability  density 
function  and  its  influence  is  then  incorporated  into  inverse  planning  through  the  use  of  a 
statistical  inference  theorem  The  technique  is  illustrated  by  using  a  hypothetical  C- 
shaped  tumor  case,  a  prostate  tumor  case  and  a  paraspinal  tumor  case  with  an  EUD-based 
model.  Considering  that  currently  available  models  for  computing  biological  effects  of 
radiation  are  simplistic,  and  the  data  they  rely  on  are  sparse  and  of  questionable  quality, 
the  proposed  technique  provides  us  with  an  effective  tool  to  minimize  the  effect  caused 
by  the  uncertainties  in  a  statistical  sense.  The  treatment  plans  so  obtained  are  generally 
less  sensitive  to  the  inter-patient  variation  and  other  types  of  uncertainties  that  may 
otherwise  influence  the  final  treatment  plan  greatly. 


METHODS  AND  MATERIALS 


Statistical  analysis-based  inverse  planning 
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The  inverse  problem  as  posed  for  IMRT  consists  of  the  determination  of  the 
beamlet  weight  vector  w  when  a  desired  plan  is  prescribed.  In  a  vectorial  form,  the  dose 
to  the  points  in  the  treatment  region  depends  upon  the  beamlet  weights  w  as 

D^=dw,  (1) 

where  d  represents  the  dose  deposition  coefficients  matrix,  expressing  the  dose  deposited 
to  any  patient  point  when  irradiated  with  a  unit  weight  beamlet.  The  total  number  of 
physically  realizable  dose  distributions  Dc  in  IMRT  is  enormous  and  increases 
exponentially  with  the  number  of  beamlets.  Inverse  planning  is  essentially  a  plan 
selection  process  from  the  vast  pool  of  physically  realizable  solutions.  In  a  recent  paper, 
Xing  et  al  introduced  a  statistical  analysis-based  inverse  planning  technique.  In  this 
approach  the  commonly  used  objective  function  is  reformulated  into  a  probability  density 
function  whose  value  gives  the  figure  of  merit  of  a  dose  distribution.  A  virtue  of  the 
approach  is  that  it  allows  us  to  obtain  solution  in  the  presence  of  uncertainties  of  the 
prescription  parameters  or  other  model  parameters  using  a  statistical  inference  technique. 
Application  of  the  technique  to  deal  with  a  system  with  a  set  of  variable  dose 
prescriptions  has  been  described  in  another  work  of  our  group  Here  we  use  the 
formalism  for  biological  modeling  based-  inverse  planning  in  the  presence  of  model 
parameter  uncertainties.  To  be  specific,  we  use  an  equivalent  uniform  dose  (EUD)-based 
objective  function  employed  by  Wu  et  al  and  discuss  the  consequences  of  the 
variation  of  model  parameter  a  and  how  to  incorporate  the  fluctuations  into  inverse 
planning  dose  optimization  to  obtain  statistically  optimal  solutions. 


5 


EUD  model  and  EUD-based  objective  function 

The  concept  of  equivalent  uniform  dose  (EUD)  for  tumor  was  originally 
introduced  by  Niemierko  as  the  biologically  equivalent  dose  that,  if  given  uniformly, 
would  lead  to  the  same  cell  kill  in  the  tumor  volume  as  the  actual  nonumiform  dose 
distribution.  Recently,  Niemierko  et  al  suggested  a  phenomenological  form 

f  \a 

(2) 

for  both  tumor  and  normal  tissues,  where  N  is  the  number  of  voxels  in  the  structure,  Di  is 
the  dose  delivered  to  the  fth  voxel,  a  is  the  tumor  or  normal  tissue-specific  parameter  that 
describes  the  dose-volume  effect.  EUD  described  in  Eq.  (2)  is  the  general  mean  of  the 
non-uniform  dose  distribution.  According  to  the  mathematic  properties  of  the  function 
for  a  =  00,  the  EUD  is  equal  to  the  maximum  dose,  and  for  a  =  -oo,  the  EUD  is  equal  to 
the  minimum  dose.  Timiors  generally  have  targe  negative  values  of  a,  whereas  serial 
critical  structures  (e.g.  spinal  cord  and  rectum)  have  large  positive  values  and  parallel 
critical  structures  that  exhibit  a  large  dose-volume  effect  (e.g.  liver,  parotids,  and  lungs) 
have  small  positive  values. 

The  objective  function  or  figure  of  merit  used  to  measure  the  goodness  of  a  dose 
distribution  and  guide  the  optimization  In  the  present  paper,  the  system  objective 
function  is  given  by 

(3) 

J 

where  the  component  subcore^  may  be  either 
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for  tumors,  or 


foAR 


^  EUD^” 


1  + 


\EUDJ 


(5) 


for  normal  tissues  and  organs  at  risk  (OARs).  EUDo  is  the  desired  dose  parameter  for  the 
target  volume  and  the  maximum  tolerable  uniform  dose  for  normal  structures.  Parameter 
n  is  akin  to  the  structure  specific  importance  factor  in  the  conventional  inverse 
planning  formalism  that  parameterizes  our  tradeoff  strategy  of  different  structure.  The 
large  n  indicates  high  importance. 


Incorporation  of  the  variation  distribution  of  the  model  parameter  into  inverse 
planning 

We  assume  that  in  the  EUD  model  varies  according  to  a  simple  Gaussian 
distribution 

=  (6) 

where  oq  is  the  mean  value,  P’^is  a  normalization  constant  and  a*  is  one  of  the  sampling 

values  of  a.  For  a  given  distribution,  the  EUD  and  the  corresponding  figure  of  merit  of  an 
IMRT  plan  vary  with  the  sampling  of  a.  We  thus  rewrite  Eqs.  (4)  and  (5)  as  conditional 
probabilities  for  a  sampled  a*: 


P^(,EUD\ak)  = 


1 


n  5 


(7) 


(eudA 
V  EUD  J 


1  + 


F,^(£[/Dlak)  = 


1 


[EUDo 


(8) 


The  objective  function  for  a  structure  m  in  the  presence  of  uncertainty  in  a  is  expressed  as 
the  summation  of  a  series  of  joint  probabilities 


PSEUD)  =  Ep  (££//)  I  aO-z’.Co.) , 

k 


(9) 


and  the  overall  objective  function  F  of  the  system  is  a  product  of  Pm(EUD)  defined  in  Eq. 
(9).  That  is 

F  =  ln(l  /P)  =  -}nYlf,(,EUD)=-'^hi'£PSEVD  |  ayP^a,) .  (10) 

m  m  k 


Optimization  method 

As  described  above,  the  uncertainties  of  model  parameters,  {ok},  are  described  by 
probability  density  functions  and  they  are  incorporated  into  the  overall  objective  function 
of  the  system  through  the  joint  probability  given  by  Eq.  (9).  To  obtain  the  optimal 
solution  in  the  presence  of  model  parameters,  all  we  need  to  do  is  to  minimize  the  overall 
objective  function  given  by  Eq.  (10). 

The  calculation  process  is  schematically  shown  in  Fig.  1.  For  the  computational 
purpose,  the  probability  density  function  for  each  structure  is  discreted  into  seven  equally 
spaced  points.  We  use  the  Fletcher-Reeves  conjugate  gradient  optimization  algorithm 
to  optimize  the  system.  But  any  other  iterative  or  stochastic  optimization  can  be  also 
employed  to  optimize  the  system.  A  common  step  in  all  optimization  algorithms  is  the 
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evaluation  of  the  objective  function  for  a  trial  beam  profiles  (or  computed  dose 
distribution),  which  is  somewhat  tedious  here  because  of  the  appearance  of  multiple  OkS 
of  the  involved  structures.  Briefly,  for  a  given  trial  beam  profiles  or  dose  distribution,  the 
evaluation  of  the  objective  function  consists  of  four  steps:  (i)  for  a  structure  m,  calculate 
the  EUD  corresponding  to  each  possible  a*;  (ii)  calculate  the  conditional  probability  for 
the  target  and  OAR  using  Eq.  (7)  and  (8),  respectively;  (iii)  sum  over  all  possible  a*  to 
obtain  the  joint  probability,  given  by  Eq.  (9);  and  (iv)  sum  over  all  structures  to  obtain  the 
overall  objective  function  value.  After  the  dose  optimization,  a  set  of  optimal  beam 
profiles  and  the  corresponding  dose  distribution  and  other  plan  indices  are  provided  for 
the  planner  to  assess  the  clinical  relevance  of  the  obtained  treatment  plan. 

Test  cases 

The  new  algorithm  was  tested  using  a  hypothetical  phantom  case  with  a  C-shaped 
target  and  two  clinical  cases  (a  prostate  case  and  a  paraspinal  tumor  boost  treatment).  The 
size  of  the  pencil  beam  defined  at  the  isocenter  was  0.5  cm.  The  configuration  of  the  C- 
shaped  tumor  case  is  shown  in  Fig.  2A.  Nine  6MV  equi-spaced  beams  were  used  for  the 
treatment  (0°,  40°,  80°,  120°,  160°,  200°,  240°,  280°,  and  320°  -  respecting  the  lEC 
convention).  The  values  of  n  and  a  in  the  EUD-based  objective  function  are  listed  in 
Table  1.  The  parameter  a  in  EUD  model  characterizes  the  dose-volume  effect  but  its 
value  is  generally  not  known  accurately  even  for  clinically  well  studied  organs.  The 
influence  of  the  uncertainty  in  the  a  value  of  a  target  or  sensitive  structure  to  the  final 
treatment  plan  was  studied  and  analyzed. 

Similar  study  was  carried  out  for  the  two  clinical  cases.  The  six  6MV  beam  angles 
used  for  the  IMRT  prostate  treatment  were  0°,  55°,  135°,  180°,  225°  and  305°.  Table  2 
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lists  some  relevant  parameters  used  for  planning  the  case.  For  the  IMRT  paraspinal  boost 
treatment,  five  6  MV  non-equally  spaced  coplanar  beams  were  placed  at  the  following 
angular  positions:  95°,  140°,  175°,  225°  and  275°.  The  target  boost  dose  was  prescribed 
to  16  Gy.  Relevant  parameters  are  listed  in  Table  2.  The  planning  goal  was  to  find  a  dose 
distribution  that  covered  the  tumor  volume  as  uniformly  as  possible,  while  maximally 
sparing  the  spinal  cord,  liver,  and  kidney. 


RESULTS  AND  DISCUSSIONS 

The  C-shaped  tumor  case 

We  first  investigated  the  behavior  of  the  system  when  the  parameter  a  of  the 
target  EUD  takes  four  different  distributions,  as  depicted  in  the  bar  charts  shown  on  the 
right  of  Fig.  3,  while  keeping  the  parameter  a  of  the  OAR  at  a  constant  ao  =  6.0  In  the 
case  shown  in  Fig.  3  dl,  the  parameter  a  takes  only  a  single  value,  ao  =  -10,  which  is  a 
simple  case  studied  by  Wu  et  al  The  optimal  plans  for  the  four  distributions  of 
parameter  a  differ  significantly,  as  indicated  by  the  target  and  OAR  DVHs  shown  in  Fig. 
3  A  and  B.  The  isodose  plot  corresponding  to  the  o-distribution  shown  in  Fig.  3  d2  is 
plotted  in  Fig.  2B. 

To  estimate  the  degree  of  sensitivity  of  the  solutions  against  a  variation  in  a,  we 
computed  the  target  EUD  and  the  objective  function,^,  as  a  function  of  parameter  a  for 
the  four  optimal  dose  distributions  under  different  types  of  uncertainty  distributions.  The 
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results  are  plotted  in  Fig.  4.  For  plan  A,  the  EUD  changes  from  65  to  71  Gy  when  a  is 
varied  from  —10  to  -70  and  to  79  Gy  when  a  is  equal  to  140.  The  objective  function  varies 
from  0.11  to  0.85  in  the  range  of  variation  in  a.  For  plan  D,  the  EUD  is  narrowed  to  a 
range  between  70  Gy  and  79  Gy.  The  EUD  variations  of  plans  B  and  C  are  similarly 
reduced.  These  results  suggest  that  the  EUD  becomes  much  less  sensitive  to  the  variation 
in  parameter  a  in  the  plans  obtained  with  some  “built-in”  distributions  in  parameter  a 
(i.e.,  plans  corresponding  to  Figs.  3  d2  to  d4). 

The  uncertainty  of  parameter  a  of  the  OAR  can  be  similarly  included  in  the  dose 
optimization  process  when  its  distribution  is  known.  In  the  second  study,  we  fixed  the 
target  EUD  parameter  a=  -10  and  allowed  the  parameter  a  of  the  OAR  to  take  four 
different  distributions  as  plotted  in  the  right  of  Fig.  5.  The  target  and  OAR  DVHs  for  the 
four  possible  scenarios  are  shown  in  A  and  B.  Once  again,  we  found  that  the  final 
solution  strongly  depends  on  the  distributions  of  the  parameter  a. 

The  maximum  doses  of  the  OAR  of  the  four  plans  vary  from  24  Gy  to  30  Gy. 
Note  that  the  doses  to  the  OAR  in  plans  B,  C  and  D  are  less  than  that  of  plan  A,  where  the 
parameter  a  is  restricted  to  a  single  value,  ao  =  6.  This  is  explainable  since  the  parameters 
a  in  plans  d2,  d3  and  d4  are  shifted  up  to  higher  values.  As  a  increases,  the  EUD  puts 
more  emphasis  on  the  high  dose  (recall  that  EUD  becomes  the  maximum  dose  when  a  = 
oo).  As  a  consequence  of  the  increased  “effective”  a  value  in  the  distributions  shown  in 
Figs.  5  d2,  d3  and  d4,  the  OAR  dose  is  improved  in  comparison  with  the  plan  obtained 
under  the  assumption  of  a  fixed  a  value  (Fig.  5  dl).  Interestingly,  the  target  DVHs  shows 
that  four  distinct  plans  have  very  similar  target  coverage.  It  is  well  known  that  in  dose 
optimization  there  is  generally  no  net  gain:  an  improvement  in  the  dose  to  a  structure  is 
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often  accompanied  by  a  dosimetrically  adverse  effect(s)  at  other  points  in  the  same  or 
different  structures.  The  result  here  suggests  that,  from  a  clinical  point  of  view,  it  is 
possible  to  have  a  great  gain  in  one  structure  with  a  little  sacrifice  in  another  structure. 
How  to  find  the  truly  optimal  tradeoff  represents  a  practical  subject  that  is  worth  of 
studying  in  the  future. 

As  can  be  expected  from  the  discussion  in  previous  paragraphs,  the  solution 
obtained  with  ao=  5  (Fig.  5  dl)  is  more  sensitive  to  a  variation  in  parameter  a.  Indeed,  as 
seen  from  Fig.  6,  the  EUD  for  this  plan  varies  from  1  Gy  to  30  Gy  when  a  is  changed 
from  -80  to  140.  On  the  other  hand,  the  EUD  changes  for  the  rest  three  situations  are 
much  less  for  the  same  variation  in  a.  The  upper  bound  of  the  EUD  is  reduced  to  26  Gy 
for  plan  d4,  24  Gy  for  plan  d2,  and  23  Gy  for  plan  d3.  The  objective  functions  of  four 
plans  show  a  similar  trend. 

The  prostate  tumor  case 

Four  IMRT  plans  with  different  types  of  pre-assumed  uncertainties  were 
generated  for  a  prostate  tumor  case  (Fig.  7A).  These  include:  (i)  The  a-parameters  for 
both  prostate  target  and  OARs  are  restricted  to  single  values  as  listed  in  Table  2.  This 
plan  serves  as  a  reference  whose  DVHs  are  shown  in  Figs.  8A-8C  as  dotted  curves;  (ii) 
Only  the  a-parameter  of  the  prostate  target  takes  a  range  of  values,  as  depicted  in  the 
right  of  Fig.  8 A;  (iii)  Only  the  o-parameter  of  the  rectum  takes  a  range  of  values,  as 
depicted  in  the  right  of  Fig.  8B;  and  (iv)  The  <3-parameters  of  both  prostate  target  and  the 
rectum  were  allowed  to  take  a  range  of  values,  as  depicted  in  the  right  of  Fig.  8C. 
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DVHs  for  the  plan  using  parameters  defined  in  Table  2  are  plotted  with  dotted 
curves  and  plans  with  the  inclusion  of  parameter  uncertainty  are  dravm  with  solid  curves 
(Fig.  8).  When  the  parameter  a  in  target  EUD  takes  a  Poisson  distribution  as  shown  in  the 
bar  chart  of  Fig.  8 A,  prostate  dose  homogeneity  is  significantly  improved.  The  minimum 
dose  increases  from  55  Gy  to  67  Gy,  and  the  maxim  dose  decreases  slightly  from  82  Gy 
to  80Gy.  However  the  volumes  receiving  radiation  dose  for  rectum,  bladder  and  normal 
tissue  all  increase  significantly  though  the  maximum  dose  remains  similar.  The 
improvement  of  the  target  coverage  and  compromise  of  OAR  sparing  is  a  natural 
outcome  of  the  competitive  requirements  for  targets  and  OARs  imposed  on  the  system. 
The  corresponding  dose  distribution  with  the  target  parameter  defined  in  the  bar  chart  A 
is  shown  in  Fig.  7B. 

Next  we  considered  the  inclusion  of  parameter  a  uncertainty  in  EUD  calculation 
in  one  of  the  critical  structures-rectum  (Fig.  8  B).  The  irradiated  rectum  volume  for  a 
dose  below  60  Gy  is  less  than  that  of  a  conventional  plan  with  the  parameter  a  fixed  at 
24.  DVHs  for  the  bladder,  normal  tissue  and  prostate  do  not  change  significantly 
compared  to  the  plan  without  inclusion  of  parameter  uncertainty. 

Lastly,  we  simultaneously  replaced  target  and  rectum  parameters  with  the 
distributions  shown  in  Fig.  8  C.  Similar  to  that  corresponds  the  prescription  of  Fig.  8  A, 
the  prostate  coverage  is  improved.  However,  the  rectum  DVH  in  this  case  is  not  worsen 
greatly  because  parameter  a  of  rectum  EUD  was  allowed  to  take  a  spectrum  of  values. 
For  bladder  and  normal  tissue,  although  their  irradiated  volumes  in  the  low  dose  region 
are  higher  than  those  of  the  conventional  plan,  the  volumes  receiving  high  doses  are 
reduced. 
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The  paraspinal  tumor  case 

Three  IMRT  plans  were  generated  for  a  paraspinal  tumor  case  (Fig.  9A).  These 
include:  (i)  The  a-parameters  for  both  prostate  target  and  OARs  are  restricted  to  single 
values  as  listed  in  Table  3.  This  plan  serves  as  a  reference  whose  DVHs  are  shown  in 
Figs.  9B-C  as  dotted  curves;  (ii)  Only  the  a-parameter  of  the  target  takes  a  range  of 
values,  as  depicted  in  the  right  of  Fig.  9B;  and  (iii)  Only  the  a-parameter  of  the  spinal 
cord  takes  a  range  of  values,  as  depicted  in  the  right  of  Fig.  9C. 

When  a  in  target  EUD  takes  a  Poisson  distribution  as  shown  in  the  bar  chart  of 
Fig.  9B,  dose  homogeneity  is  slightly  improved.  However  this  is  achieved  at  the  expense 
of  more  irradiation  to  the  cord.  The  inclusion  of  parameter  a  uncertainty  in  EUD 
calculation  in  the  spinal  cord  (Fig.  9C)  reduced  the  maximum  cord  dose  by  100  cGy.  The 
DVHs  of  the  target,  kidney,  and  liver  were  not  changed  significantly  compared  to  the 
plan  without  inclusion  of  parameter  uncertainty. 

The  influence  of  various  uncertainties  on  the  patient  treatment  has  been  a  subject 
of  intense  study.  Fenwick  and  Nahum  have  included  the  model  parameter  uncertainty 
with  a  standard  deviation  when  calculating  the  NTCP  of  rectum  Similarly,  the 
inclusion  of  uncertainties  in  the  patient  setup  and  dose  calculation  has  also  been 
demonstrated  Deasy  et  al  have  used  a  bootstrap-based  method  to  estimate  the 
influence  of  biological  parameter  uncertainties  on  predicting  long-term  salivary  function 
.  The  statistical  method  proposed  here  provides  a  general  framework  to  include  various 
uncertainties  in  the  dose  optimization  process.  With  minor  modification,  the  technique 
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can  be  extended  to  derive  statistically  optimal  solutions  in  the  presence  of  other  types  of 
uncertainties. 

As  can  be  intuitively  imagined,  the  inclusion  of  a-distribution  will  definitely 
change  the  final  dose.  Whether  it  will  improve  or  worsen  the  final  dose  distribution  will 
generally  depend  on  the  specific  form  of  the  a-distribution,  and  also  the  metric  used  to 
judge  the  goodness  of  a  plan.  If  the  original  EUD-based  objective  function  is  used  as  the 
sole  metric  for  the  judgment,  the  inclusion  of  a-distribution  may  make  the  plan  worse. 
However,  clinical  decision-making  is  not  made  by  a  single  function  and  a  “worse”  plan 
judged  by  the  EUD-objective  function  may  turn  out  to  be  clinically  more  favorable.  In 
other  words,  there  is  a  gap  between  mathematical  dose  optimization  and  clinical  decision¬ 
making.  The  study  seems  to  suggest  that,  while  it  is  generally  true  that  there  is  no  net 
gain  in  dose  optimization  it  is  important  to  develop  a  method  that  is  capable  of 
optimizing  not  only  the  objective  function  but  also  the  next  level  of  decision-making. 
This  kind  of  optimization  will  allow  us  to  find  the  solution  that  may  sacrifice  a  little  (i.e., 
clinically  insignificant)  in  one  or  a  few  structures  but  gain  a  lot  in  other  structures. 

CONCLUSIONS 

We  have  proposed  and  implemented  a  technique  for  incorporating  biological 
model  parameter  uncertainties  into  inverse  treatment  planning.  The  formalism  is  quite 
general  and  does  not  prerequisite  the  specific  form  of  uncertainty  distributions  of  the 
involved  model  parameters.  By  including  model  parameter  uncertainties,  the  final 
solution  becomes  more  robust  and  the  treatment  outcome  will  be  less  likely  influenced  by 
inter-patient  variation  of  biological  characteristics.  With  the  increasing  interest  in 
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radiation  therapy  community  to  use  biologically  based  models  for  treatment  planning, 
this  work  provides  an  effective  way  to  better  account  for  the  known  uncertainties  in  the 
model  parameters  and  allows  us  to  maximally  utilize  the  available  radiobiology 
knowledge  to  facilitate  patient  care. 
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FIGURE  CAPTIONS 

Figure  1.  The  flow  chart  of  the  optimization  process  with  the  inclusion  of  model 
parameter  uncertainty. 

Figure  2.  (A)  The  sketch  of  the  hypothetical  case  with  C-shaped  target  and  the  beam  set¬ 
up  for  dose  optimization.  (B)  The  dose  distribution  corresponding  to  the  parameters  listed 
in  Tab.l  and  the  probabilistic  distribution  shown  in  Fig.  3B. 

Figure  3.  The  target  and  OAR  DVHs  of  four  optimal  plans  when  parameter  a  is  a  fixed 
value  (bar  chart  dl)  and  varies  according  to  three  different  probabilistic  distributions  (bar 
chart  d2,  d3  and  d4). 

Figure  4.  The  EUD  of  the  target  and  objective  function  when  parameter  a  is  prescribed 
according  to  Fig.  3. 

Figure  5.  The  target  and  OAR  DVHs  of  four  optimal  plans  when  parameter  <at  is  a  fixed 
value  (bar  chart  dl)  and  varies  according  to  three  different  probabilistic  distributions  (bar 
chart  d2,  d3  and  d4). 

Figure  6.  The  EUD  of  the  OAR  and  objective  function  when  parameter  a  is  prescribed 
according  to  Fig.  5. 
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Figure  7.  A  transverse  slice  showing  the  anatomical  structures  delineated  for  the  prostate 
tumor  (A)  and  the  corresponding  optimized  dose  distribution  with  the  parameters  listed  in 
Tab.  2  and  the  probabilistic  distribution  shown  in  Fig.  8B. 

Figure  8.  DVHs  for  a  prostate  cancer  case  using  the  conventional  optimization  with  fixed 
flf-value  (dotted  line)  and  the  newly  proposed  approach  with  the  inclusion  of  model 
parameter  uncertainty  (solid  line).  (A)  Only  the  a-parameter  for  the  target  is  assigned 
with  a  probabilistic  distribution;  (B)  Only  the  a-parameter  for  the  OAR  is  assigned  with  a 
probabilistic  distribution;  (C)  Uncertainties  in  the  a-parameter  are  introduced  for  both  the 
target  and  OAR. 

Figure  9.  (A)  A  transverse  slice  showing  the  anatomical  structures  for  a  paraspinal  case; 
(B)  DVHs  when  the  a-parameter  for  the  target  is  assigned  with  a  probabilistic 
distribution;  (C)  DVHs  when  the  a-parameter  for  the  OAR  is  assigned  with  a 
probabilistic  distribution. 
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Table  1 .  The  conventional  EUD-based  optimization  parameter  for  the  hypothetical  IMRT 
treatment  of  a  C-shaped  tumor. 


PTV 

PTV* 

OAR 

NT 

a 

o 

o 

1 

io.o 

6.0 

6.0 

72 

76 

35 

35 

n 

20 

6 

6 

*  Contains  parameters  for  the  target  treated  as  virtual  normal  tissue  to  limit  dose 
inhomogeniety. 


Table  2.  The  conventional  EUD-based  optimization  parameter  for  prostate  cancer. 


PTV 

PXV* 

Bladder 

Rectum 

NT 

a 

-10.0 

10.0 

6.0 

24 

6.0 

Mmmsm 

72 

76 

35 

35 

35 

n 

20 

20 

6 

6 

6 

*  Contains  parameters  for  the  target  treated  as  virtual  normal  tissue  to  limit  dose 
inhomogeniety. 


Table  3.  The  conventional  EUD-based  optimization  parameter  for  paraspinal  tumor. 


PTV 

PTV* 

Liver 

a 

-10.0 

10.0 

6.0 

6.0 

6.0 

lOSiKSSH 

16 

17 

12 

6.4 

OO 

n 

20 

20 

6 

6 

6 

*  Contains  parameters  for  the  target  treated  as  virtual  normal  tissue  to  limit  dose 
inhomogeniety. 
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Abstract 

Radiobiological  treatment  planning  depends  not  only  on  the  accuracy  of  the  models  describing  the  dose-response  relation  of 
different  tumors  and  normal  tissues  but  also  on  the  accuracy  of  tissue  specific  radiobiological  parameters  in  these  models.  Whereas 
the  general  formalism  remains  the  same,  different  sets  of  model  parameters  lead  to  different  solutions  and  thus  critically  determine 
the  final  plan.  Here  we  describe  an  inverse  planning  formalism  with  inclusion  of  model  parameter  uncertainties.  This  is  made 
possible  by  using  a  statistical  analysis-based  fi-ameset  developed  by  our  group.  In  this  formalism,  the  uncertainties  of  model 
parameters,  such  as  the  parameter  a  that  describes  tissue-specific  effect  in  EUD  model,  are  expressed  by  probability  density 
function  and  are  included  in  the  dose  optimization  process.  We  found  that  the  final  solution  strongly  depends  on  distribution 
functions  of  the  model  parameters.  Considering  that  currently  available  models  for  computing  biological  effects  of  radiation  are 
simplistic,  and  the  clinical  data  used  to  derive  the  models  are  sparse  and  of  questionable  quality,  the  proposed  technique  provides  us 
with  an  effective  tool  to  minimize  the  effect  caused  by  the  uncertainties  in  a  statistical  sense.  With  the  incorporation  of  the 
imcertainties,  the  technique  has  potential  for  us  to  maximally  utilize  the  available  radiobiology  knowledge  for  better  IMRT 
treatment. 


Keywords 

inverse  planning,  dose  optimization,  biological  models,  IMRT 

Introduction 

Most  IMRT  optimization  systems  at  present  use  dose  and/or 
dose  volume-based  objective  functions  [1, 2],  which  guide  the 
IMRT  planning  by  imposing  a  penalty  according  to  the 
difference  between  the  computed  and  prescribed  doses.  A 
well-known  drawback  of  the  dose-based  inverse  planning  is 
that  the  nonlinear  dose  response  of  tumor  or  normal  structures 
is  not  fully  considered.  A  number  of  mathematical  models 
have  been  developed  over  the  years  to  better  describe  the 
biological  effect  of  radiation,  which  include  tumor  control 
probability  (TCP)  [3],  normal  tissue  complication  probability 
(NTCP)  [4],  equivalent  uniform  dose  (EUD)  [5]  and  the 
probability  of  uncomplicated  tumor  control  (P+)  [6].  In 
parallel  to  these  modeling  efforts,  considerable  works  have 
also  been  done  to  use  these  biological  models  to  construct 
more  meaningful  objective  functions  for  therapeutic  dose 
optimization  [7] . 

Generally  speaking,  radiobiological  formalism  involves  the 
use  of  model  parameters  that  are  of  considerable  uncertainty 
[8].  For  instance,  the  radiosensitivity  a  of  Webb’s  TCP  model 
varies  from  0.157  Gy-1  to  0.090  Gy-1  when  model 
parameters  were  fit  to  103  patients’  data  [3].  Biological 
‘margins’  have  been  used  to  account  for  the  variability  in 
radiation  sensitivity.  This  method  assumes  the  patient  to  be 
more  sensitive  than  the  mean  value  for  normal  tissues  and 
more  resistant  for  the  tumor  [9].  K&ver  et  al  proposed  a 
stochastic  optimization  to  account  for  clinical  uncertainties, 
including  the  vaiying  radiosensitivity  [10,  1 1].  The  objective 
function  was  constructed  based  on  a  linear  quadratic  Poisson 


model  which  approximates  the  probability  of  curing  the 
patient  or  inflicting  injury.  Two  parameters  in  the  model 
could  be  calculated  if  the  standard  deviation  of  dose  per 
faction  was  known.  The  optimization  was  thus  executed 
corresponding  to  different  standard  deviations. 

We  have  recently  presented  a  general  statistical  analysis- 
based  inverse  planning  fi'amework  [12,  13]  and  applied  it  to 
investigate  the  influence  of  model  parameter  imcertainties  in 
biologically  based  dose  optimization  [14].  The  purpose  of 
this  paper  is  to  provide  a  detailed  description  of  the  technique 
and  addresses  several  important  issues  related  to  the  dose 
optimization  in  the  presence  of  model  parameter 
uncertainties.  In  our  approach,  the  uncertainty  of  a  model 
parameter  is  quantified  by  a  probability  density  function  and 
its  influence  is  then  incorporated  into  inverse  planning 
through  the  use  of  a  statistical  inference  theorem  [12].  The 
technique  is  illustrated  by  using  a  hypothetical  C-shaped 
tumor  case  and  a  prostate  tumor  case  with  an  EUD-based 
model.  Considering  that  currently  available  models  for 
computing  biological  effects  of  radiation  are  simplistic,  and 
the  data  they  rely  on  are  sparse  and  of  questionable  quality, 
the  proposed  technique  provides  us  with  an  effective  tool  to 
minimize  the  effect  caused  by  the  uncertainties  in  a  statistical 
sense.  The  treatment  plans  so  obtained  are  generally  less 
sensitive  to  the  inter-patient  variation  and  other  types  of 
uncertainties  that  may  otherwise  influence  the  final  treatment 
plan  greatly. 


^  Material  and  Methods 


Statistical  analysis-based  inverse  planning 


The  inverse  problem  as  posed  for  IMRT  consists  of  the 
determination  of  the  beamlet  weight  vector  w  when  a 
desired  plan  is  prescribed.  In  a  vectorial  form,  the  dose  to 
the  points  in  the  treatment  region  depend  upon  the  beamlet 
weights  w  as 

D^=dw,  (1) 

where  d  represents  the  dose  deposition  coefiRcients  matrix, 
expressing  the  dose  deposited  to  any  patient  point  when 
irradiated  with  a  unit  weight  beamlet.  The  total  number  of 
physically  realizable  dose  distributions  Dc  in  IMRT  is 
enormous  and  increases  exponentially  with  the  number  of 
beamlets.  Inverse  planning  is  essentially  a  plan  selection 
process  from  the  vast  pool  of  physically  realizable 
solutions.  In  a  recent  paper,  Xing  et  al  [12]  introduced  a 
statistical  analysis-based  inverse  planning  technique.  In  this 
approach  the  commonly  used  objective  fimction  is 
reformulated  into  a  probability  density  function  whose 
value  gives  the  figure  of  merit  of  a  dose  distribution.  A 
virtue  of  the  approach  is  that  it  allows  us  to  obtain  solution 
in  the  presence  of  uncertainties  of  the  prescription 
parameters  or  other  model  parameters  using  a  statistical 
inference  technique.  Application  of  the  technique  to  deal 
with  a  system  with  a  set  of  variable  dose  prescriptions  has 
been  described  in  another  work  of  our  group  [13].  Here  we 
use  the  formalism  for  biological  modeling  based-  inverse 
planning  in  the  presence  of  model  parameter  uncertainties. 
To  be  specific,  we  use  an  equivalent  uniform  dose  (EUD)- 
based  objective  fimction  employed  by  Wu  et  al  [7]  and 
discuss  the  consequences  of  the  variation  of  model 
parameter  a  and  how  to  incorporate  the  fluctuations  into 
inverse  planning  dose  optimization  to  obtain  statistically 
optimal  solutions. 


EUD  model  and  EUD-based  objective  function 
The  concept  of  equivalent  uniform  dose  (EUD)  for  tumor 
was  originally  introduced  by  Niemierko  as  the  biologically 
equivalent  dose  that,  if  given  uniformly,  would  lead  to  the 
same  cell  kill  in  the  tumor  volume  as  the  actual 
nonumiform  dose  distribution.  Recently,  Niemierko  et  al 
suggested  a  phenomenological  form  [5]: 
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(2) 


for  both  tumor  and  normal  tissues,  where  N  is  the  number 
of  voxels  in  the  structure,  a  is  the  tumor  or  normal  tissue- 
specific  parameter  that  describes  the  dose-volume  effect. 
EUD  described  in  Eq.  (2)  is  the  general  mean  of  the  non- 
uniform  dose  distribution.  According  to  the  mathematic 
properties  of  the  function  [15],  for  a  =  oo,  the  EUD  is  equal 
to  the  maximal  dose,  and  for  a  =  -c»,  the  EUD  is  equal  to 
the  minimum  dose.  Tumors  generally  have  large  negative 
values  of  a,  whereas  serial  critical  structures  (e.g.  spinal 
cord  and  rectum)  have  large  positive  values  and  parallel 
critical  structures  that  exhibit  a  large  dose-volume  effect 
(e.g.  liver,  parotids,  and  lungs)  have  small  positive  values. 


The  objective  fimction  or  figure  of  merit  used  to 
measure  the  goodness  of  a  dose  distribution  or  the 
corresponding  EUD  is  given  by  [7] 
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where  the  component  subcore^  may  be  either 
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for  normal  tissues  and  organs  at  risk  (OARs).  EUDo  is  the 
desired  dose  parameter  for  the  target  volume  and  the 
maximal  tolerable  uniform  dose  for  normal  structures. 
Parameter  n  is  akin  to  the  structure  specific  importance 
factor  [16]  in  the  conventional  inverse  planning  formalism 
that  parameterizes  our  tradeoff  strategy  of  different 
structure. 


Incorporation  of  the  variation  distribution  of  the  model 
parameter  into  inverse  planning 

We  assume  that  a*  in  the  EUD  model  varies  according  to  a 
simple  Gaussian  distribution 

=  (6) 


where  ao  is  the  mean  value,  P\  is  a  normalization  constant 

and  o*  is  one  of  the  sampling  values  of  a.  For  a  given 
distribution,  the  EUD  and  the  corresponding  figure  of  merit 
of  an  IMRT  plan  vary  with  the  sampling  of  a.  We  thus 
rewrite  Eqs.  (4)  and  (5)  as  condition^  probabilities  for  a 
sampled  a*: 
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The  objective  function  for  a  structure  m  in  the  presence  of 
uncertainty  in  a  is  expressed  as  the  summation  of  a  series  of 
joint  probabilities 


k 


and  the  overall  objective  function  P  of  the  system  is  a 
product  of  PJiEUD)  defined  in  Eq.  (9).  That  is 


F=]n(l/P)=-hn7i:(^)=-IlnX7i:.(£t^l^^^K(«t) 
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Results  and  Discussion 


The  C-shaped  tumor  case 


-  We  illustrate  the  proposed  technique  by  planning  a 
hypothetical  IMRT  phantom  case  with  a  C-shaped  target 
next  to  a  circular  critical  structure.  Nine  6MV  equi-spaced 
beams  were  used  for  the  treatment  (0°,  40°,  80°,  120°,  160°, 
200°,  240°,  280°,  and  320°  -  respecting  the  lEC 
convention). 

The  parameter  a  in  EUD  model  characterizes  the  dose- 
volume  effect  but  its  value  is  generally  not  known 
accurately  even  for  clinically  well  studied  organs.  We  first 
investigated  the  behavior  of  the  system  when  the  parameter 
a  of  the  target  EUD  takes  four  different  distributions,  as 
depicted  in  the  bar  charts  shown  on  the  right  of  Fig.  1, 
while  keeping  the  parameter  a  of  the  OAR  at  a. constant  ao= 
6.0.  In  the  case  shown  in  Fig.  lA,  the  parameter  a  takes 
only  a  single  value,  ao  = -10,  which  is  a  simple  case  studied 
by  Wu  et  al  [7].  The  rest  of  Fig.  1  shows  three 
representative  types  of  distributions  of  the  EUD  parameter 
a.  For  each  of  these  situations  we  carried  out  the  dose 
optimization  calculation  using  the  method  outlined  in  the 
last  section. 

The  optimal  plans  for  the  four  distributions  of  par^eter  a 
differ  significantly,  as  indicated  by  the  target  and  OAR 
DVHs  shown  in  Fig.  2.  To  estimate  the  degree  of  sensitivity 
of  the  solutions  against  a  variation  in  a,  we  computed  the 
target  EUD  and  the  objective  function,;^,  as  a  function  of 
parameter  a  for  the  four  optimal  dose  distributions  under 
different  types  of  uncertainty  distributions.  The  results  are 
plotted  in  the  left  panel  of  Fig.  1.  For  plan  A,  the  EUD 
changes  from  65  to  71  Gy  when  a  is  varied  fi'om  -10  to  -70 
and  to  79  Gy  when  a  is  equal  to  140.  The  objective  function 
varies  from  0.11  to  0.85  in  the  range  of  variation  in  a.  For 
plan  D,  the  EUD  is  narrowed  to  a  range  between  70  Gy  and 
79  Gy.  The  EUD  variations  of  plans  B  and  C  are  similarly 
reduced.  These  results  suggest  that  the  EUD  becomes  much 
less  sensitive  to  the  variation  in  parameter  a  in  the  plans 
obtained  with  some  “built-in”  distributions  in  parameter  a 
(i.e.,  plans  corresponding  to  Figs.  1B-3D). 

The  uncertainty  of  parameter  a  of  the  OAR  can  be  similarly 
included  in  the  dose  optimization  process  when  its 
distribution  is  known.  In  the  second  study,  we  fixed  the 
target  EUD  parameter  a=  -10  and  allowed  the  parameter  a 
of  the  OAR  to  take  four  different  distributions  as  plotted  in 
the  right  of  Fig.  3.  Once  again,  we  found  that  the  final 
solution  strongly  depends  on  the  distributions  of  the 
parameter  a. 

In  Fig.  4  we  plot  the  target  and  OAR  DVHs  for  the  four 
possible  scenarios  shown  in  Fig.  3.  The  maximum  doses  of 
the  OAR  of  the  four  plans  vaiy  from  24  Gy  to  30  Gy.  Note 
that  the  doses  to  the  OAR  in  plans  B,  C  and  D  are  less  than 
that  of  plan  A,  where  the  parameter  a  is  restricted  to  a 
single  value,  ao-  6.  This  is  explainable  since  the  parameters 
a  in  plans  B,  C  and  D  are  shifted  up  to  higher  values.  As  a 
increases,  the  EUD  puts  more  emphasis  on  the  high  dose 
(recall  that  EUD  becomes  the  maximum  dose  when  a  =  00). 
As  a  consequence  of  the  increased  “effective”  a  value  in  the 
distributions  shown  in  Figs.  3  B,  C  and  D,  the  OAR  dose  is 
improved  in  comparison  with  the  plan  obtained  under  the 
assumption  of  a  fixed  a  value  (Fig.  3A).  Interestingly,  the 
target  DVHs  shows  that  four  distinct  plans  have  very 
similar  target  coverage.  It  is  well  known  that  in  dose 
optimization  there  is  generally  no  net  gain;  an  improvement 


in  the  dose  to  a  structure  is  often  accompanied  by  a 
dosimetrically  adverse  effect(s)  at  other  points  in  the  same 
or  different  structures.  The  result  here  suggests  that,  from  a 
clinical  point  of  view,  it  is  possible  to  have  a  great  gain  in 
one  structure  with  a  little  sacrifice  in  another  structure. 
How  to  find  the  truly  optimal  tradeoff  represents  a  practical 
subject  that  is  worth  of  studying  in  the  future. 

As  can  be  expected  from  the  discussion  in  previous 
paragraphs,  the  solution  obtained  with  ao=  6  (Fig.  3  a)  is 
more  sensitive  to  a  variation  in  parameter  a.  Indeed,  as  seen 
from  Fig.  3,  the  EUD  for  this  plan  varies  from  1  Gy  to  30 
Gy  when  a  is  changed  from  -80  to  140.  On  the  other  hand, 
the  EUD  changes  for  the  rest  three  situations  are  much  less 
for  the  same  variation  in  a.  The  upper  bound  of  the  EUD  is 
reduced  to  26  Gy  for  plan  D,  24  Gy  for  plan  B,  and  23  Gy 
for  plan  C.  The  objective  functions  of  four  plans  show  a 
similar  trend. 

Conclusions 

We  have  proposed  and  implemented  a  technique  for 
incorporating  biological  model  parameter  imcertainties  into 
inverse  treatment  planning.  The  formalism  is  quite  general 
and  does  not  prerequisite  the  specific  form  of  uncertainty 
distributions  of  the  involved  model  parameters.  By 
including  model  parameter  uncertainties,  the  final  solution 
becomes  more  robust  and  the  treatment  outcome  will  be 
less  likely  influenced  by  inter-patient  variation  of  biological 
characteristics.  With  the  increasing  interest  in  radiation 
therapy  community  to  use  biologically  based  models  for 
treatment  planning,  this  work  provides  ah  effective  way  to 
better  account  for  file  known  uncertainties  in  the  model 
parameters  and  allows  us  to  maximally  utilize  the  available 
radiobiology  knowledge  to  facilitate  patient  care. 
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Abstract 

With  the  development  of  functional  imaging  techniques  and  intensity  modulated  radiation 
therapy  (IMRT),  there  is  growing  interest  in  image-guided  IMRT  to  produce  customized  3D 
dose  distributions  in  accordance  with  the  patient  specific  biological  requirements.  MRSI  is 
one  of  the  most  promising  functional  imaging  modalities  and  has  been  applied  to  brain, 
breast,  and  prostate  cancer  imaging.  In  practice,  however,  the  MRSI  data  do  not  always 
accurately  reflect  the  actual  metabolic  level  over  the  entire  imaging  volume  due  to  our  limited 
ability  to  shim  near  air-filled  cavities  and/or  the  strong  dependence  of  the  surface  coil  SNR  on  the 
spatial  position.  In  this  work,  we  provide  an  algorithm  to  numerically  incorporate  the  spectral 
uncertainties  into  MRSI-guided  IMRT  treatment  planning.  Assuming  that  the  fluctuation  of  the 
spectral  activity  or  the  prescribed  dose  EUDo(n)  at  voxel  n  is  specified  by  a  probability 
chstribution  P„(EUDo),  we  treat  P„(EUDo)  as  a  priori  variable  distribution  and  constmct  an 
objective  function  based  on  the  statistical  inference  technique.  The  algorithm  is  used  to  plan  a 
phantom  case  with  hypothetical  functional  distributions  and  a  brain  tumor  treatment  with 
incorporation  of  MRSI  and  the  corresponding  confidence  map.  The  results  indicated  that  the 
proposed  technique  is  capable  of  producing  deliberately  non-uniform  dose  distributions 
consistent  with  the  MRSI  and  its  spatial  uncertainty  distribution.  Considering  that  currently 
available  functional  image  data  are  not  completely  reliable  and  that  missing  spectral  data 
may  occur  frequently,  the  proposed  technique  provides  us  with  an  effective  tool  to  minimize 
the  effect  and  generate  statistically  optimal  treatment  plans. 


Including  Metabolic  Uncertainty  into  Proton  MR  Spectroscopic  Imaging 
(MRSI)-Guided  Inverse  Treatment  Planning 


Introduction 

MRSI  is  one  of  the  most  important  functional  imaging  modalities  and  has  been  applied  to  brain,  breast,  and 
prostate  cancer  imaging.  The  modality  can  not  only  be  used  to  more  accurately  delineate  the  tumor  target  but  also 
reveals  tumor  biology  distribution  that  allows  us  to  identify  high/low  tumor  burden  regions.  The  imaging  data  are 
thus  potentially  useful  to  guide  radiation  therapy  treatment  planning  to  produce  deliberately  non-uniform  dose 
distribution  that  selectively  boosts  the  high  tumor  burden  regions.  In  practice,  however,  the  MRSI  data  do  not  always 
accurately  reflect  the  actual  metabolic  level  over  the  entire  imaging  volume  due  to  our  limited  ability  to  shim  near  air- 
filled  cavities  and/or  the  strong  dependence  of  the  surface  coil  SNR  on  the  spatial  position.  In  this  work,  we  provide  an 
algorithm  to  numerically  incorporate  the  spectral  uncertainties  (confidence  map)  into  MRSI-guided  IMRT  treatment 
planning.  The  new  formalism  is  based  on  the  Bayesian  statistical  inference  theorem.  For  illustration  purpose,  a  EUD 
model  is  used  as  the  plan  ranking  function.  The  method  is  applied  to  study  a  phantom  case  with  a  few  hypothetical 
functional  distributions  and  a  brain  tumor  case  with  inhomogeneous  tissue  abnormal  levels  indicated  by  a  ratio  of 
Choline  and  NAA  firom  MRSI. 

Methods 

The  equivalent  uniform  dose  (EUD)  is  the  biologically  equivalent  dose  that,  if  given  uniformly,  would  lead 
to  the  same  cell  kill  in  the  structure  as  the  actual  non-uniform  dose  distribution.  It  is  defined  as: 

for  both  tumor  and  normal  tissues,  where  A  is  the  number  of  voxels  in  a  structure,  a  is  the  tissue-specific  parameter 
that  describes  the  dose- volume  effect.  Assuming  EUDq  (n)  is  the  voxel  desired  dose  parameter  and  it  is  linearly 
related  to  tissue  metabolic  level  M(n)  (Cho/NAA  in  brain  tumor).  We  postulate 

EUD^(n)  =  EUD;(n)  +  kM(n) 

for  a  tumor,  where  EUDq  (n)  is  the  desired  dose  at  voxel  n,  EUDQ‘'(n)  is  the  conventional  prescription  dose  and  k  is 

an  empirical  coefficient.  Similarly  EUDq  (n)  for  a  critical  structure  is  defined  linearly  proportional  to  a  functional 
importance  factor.  In  this  way,  the  uncertainty  of  measured  biological  data  can  be  projected  into  a  EUDq  distribution 
function.  Given  a  value  of  EUDq,  our  preference  over  the  occurrence  of  the  EUD  can  be  expressed  as  a  conditional 
probability, 

P^iEUD\EUDj - 
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for  target  and  organ  at  risk,  respectively.  The  uncertainty  in  ElJDo/M(n)  can  be  cast  into  the  objective  function  or  the 
preference  function  of  the  system  based  on  Bayesian  theorem.  The  preference  can  be  modeled  as  the  sunnnation  of  a 
series  of  joint  probabilities: 

P  (EUD)  =  ^  P  (EUD  I  EUD^^  ).P  (PC/D  ^ ) 
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The  overall  preference  function  PofN  voxel  system  is  a  product  of  Pn(EUD). 

Results 

The  method  was  applied  to  an  IMRT  treatment  of  malignant  glioma.  The  target  metabolic  map  was 
discretized  into  three  levels  ^ig.  A)  based  on  the  values  of  choline/NAA  ratio.  In  a  conventional  plan,  the  inner 
region  with  highest  abnormality  (target  1)  was  prescribed  64  Gy,  the  middle  region  (target  2)  54Gy  and  the  external 
region  (target  3)  44  Gy.  Because  the  measurement  of  choline/NAA  has  great  uncertainty,  we  replaced  the  fixed 


**  i  ; 


value  desired  dose  (EUDo)  with  a  probability  density  based  dose  prescription.  For  instance,  the  target  1  was 
prescribed  58, 60, 62, 64, 66, 68, and  70Gy  with  probabilities  12%,  13%, 16%, 18%,  16%,  13%  and  12%  respectively.  A 
dose  distribution  is  shown  in  Fig.  B  when  three  targets  are  all  prescribed  with  a  Gaussian  type  distribution. 

The  objective  functions  are  lower  than  those  of  the  conventional  plan  (Fig.  C).  The  corresponding  DVH  shows  this 
type  of  distribution  could  result  in  slight  under  dose  (Fig.  D),  suggesting  a  higher  dose  prescription  may  be  needed 
when  measurement  uncertainty  is  of  a  Gaussian  form.  Target  DVHs  highly  depend  on  the  distribution  of  the 
abnormity  levels.  Assuming  target  3  needs  to  be  prescribed  with  higher  probabilities  for  doses  less  than  44  Gy,  and 
target  1  needs  to  be  prescribed  with  higher  probabilities  for  doses  over  64  Gy,  we  found  the  resultant  DVHs  were 
significantly  different  than  those  of  a  conventional  plan  (Fig.  E), 

Conclusions 

We  have  used  functional  MRSI  metabolic  data  to  guide  the  design  of  IMRT  treatment  plan.  The  uncertainty 
represented  in  functional  imaging  has  been  integrated  into  the  dose  optimization.  Using  this  algorithm,  the  effects  of 
the  MRSI  spectral  uncertainty  can  be  minimized  in  a  statistical  sense  and  functional  data  can  be  used  more 
efficiently  and  accurately. 

This  work  was  supported  in  part  by  the  US  Army  Medical  Research  Grant. 
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A)  Three  targets  (1,2  and  3)  with  different  abnormality  indicted  by  MRSI.  B)  Dose  distribution  with  inclusion  of 
Gaussian  type  uncertainty.  Target  1,  58  to  70Gy,  target  2,  48  to  60Gy  and  target  3,  38  to  50Gy,  all  with  2  Gy  intervals. 
The  probability  distributions  are  12%,  13%,  16%,  18%,  16%,  13%  and  12%  for  seven  discretized  dose  levels.  C)  Target 
objective  functions  for  dose  prescription  stated  in  B.  D)  Target  DVHs  corresponding  to  dose  prescription  stated  in  B.  E) 
Target  DVHs  corresponding  to  probability  distributions:  target  1,1%,  2%,  3%,  39%,  21%,  18%,  and  16%;  target  2, 

12%,  13%,  16%,  18%,  16%,  13%  and  12%;  target  3,  16%,  18%,  21%,  39%,  3%,  2%,  and  1%,  Discretized  dose  levels 
are  the  same  as  in  B. 


